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Introduction to Statistical Signal Processing 


Digital Signal Processing 


e Digital = sampled, discrete-time, quantized 
e Signal = waveform, sequnce of measurements or observations 
e Processing = analyze, modify, filter, synthesize 


Examples of Digital Signals 


e sampled speech waveform 
e "pixelized" image 
¢ Dow-Jones Index 


DSP Applications 


e Filtering (noise reduction) 
e Pattern recognition (speech, faces, fingerprints) 
e Compression 


A Major Difficulty 


In many (perhaps most) DSP applications we don't have complete or perfect 
knowledge of the signals we wish to process. We are faced with many 
unknowns and uncertainties. 


Examples 


e noisy Measurements 

e unknown signal parameters 

e noisy system or environmental conditions 

e natural variability in the signals encountered 


Functional Magnetic Resonance Imaging 
[missing resource: FMRI.png] 


Challenges are 
measurement noise and 
intrinsic uncertainties in 

signal behavior. 


How can we design signal processing algorithms in the face of such 
uncertainty? 


Can we model the uncertainty and incorporate this model into the design 
process? 


Statistical signal processing is the study of these questions. 


Modeling Uncertainty 


The most widely accepted and commonly used approach to modeling 
uncertainty is Probability Theory (although other alternatives exist such as 
Fuzzy Logic). 


Probability Theory models uncertainty by specifying the chance of 
observing certain signals. 


Alternatively, one can view probability as specifying the degree to which 
we believe a signal reflects the true state of nature. 


Examples of Probabilistic Models 


¢ errors in a measurement (due to an imprecise measuring device) 
modeled as realizations of a Gaussian random variable. 

¢ uncertainty in the phase of a sinusoidal signal modeled as a uniform 
random variable on [0, 27). 

¢ uncertainty in the number of photons stiking a CCD per unit time 
modeled as a Poisson random variable. 


Statistical Inference 


A Statistic is a function of observed data. 


Example: 
Suppose we observe WN scalar values 21,..., yy. The following are 
Statistics: 


er=4 S~_, &n (sample mean) 
° 21,-..,Xy (the data itself) 

e min {z,,..., 2} (order statistic) 
° (x17 — £2 sin(x3), e '"!*)) 


A Statistic cannot depend on unknown parameters. 


Probability is used to model uncertainty. 


Statistics are used to draw conclusions about probability models. 


probability 


Population 
of measured 
“random” signal 
signals 


statistics 


Probability models our uncertainty about signals we may observe. 


Statistics reasons from the measured signal to the population of possible 
signals. 


Statistical Signal Processing 
e Step 1 Postulate a probability model (or models) that reasonably 
capture the uncertainties at hand 
e Step 2 Collect data 


e Step 3 Formulate statistics that allow us to interpret or understand our 
probability model(s) 


In this class 


The two major kinds of problems that we will study are detection and 
estimation. Most SSP problems fall under one of these two headings. 


Detection Theory 


Given two (or more) probability models, which on best explains the signal? 


Examples 


1. Decode wireless comm signal into string of 0's and 1's 
2. Pattern recognition 


© voice recognition 
o face recognition 
o handwritten character recognition 


3. Anomaly detection 


© radar, sonar 
© irregular, heartbeat 
© gamma-ray burst in deep space 


Estimation Theory 


If our probability model has free parameters, what are the best parameter 
settings to describe the signal we've observed? 


Examples 


1. Noise reduction 
2. Determine parameters of a sinusoid (phase, amplitude, frequency) 
3. Adaptive filtering 


o track trajectories of space-craft 
© automatic control systems 


o channel equalization 


4. Determine location of a submarine (sonar) 
5. Seismology: estimate depth below ground of an oil deposit 


Example: 


Detection Example 


Suppose we observe WN tosses of an unfair coin. We would like to decide 
which side the coin favors, heads or tails. 


e Step 1 Assume each toss is a realization of a Bernoulli random 
variable. 


Pr[Heads] = p = 1 — Pr[Tails| 
Must decide p = 4 Sa — a 
e Step 2 Collect data 71,...,2y 
x; = 1 = Heads 
x; —0= Tails 


e Step 3 Formulate a useful statistic 


N 
by er 


= 


lings Be 2, guess p = +.1fk > + guess p = 4. 


Example: 
Estimation Example 


Suppose we take N measurements of a DC voltage A with a noisy 
voltmeter. We would like to estimate A. 


e Step 1 Assume a Gaussian noise model 


L, =At+wu, 


where w, ~ (0,1). 
e Step 2 Gather data 71,...,2y 
e Step 3 Compute the sample mean, 


and use this as an estimate. 


In these examples (({link] and [link]), we solved detection and estimation 
problems using intuition and heuristics (in Step 3). 


This course will focus on developing principled and mathematically 
rigorous approaches to detection and estimation, using the theoretical 
framework of probability and statistics. 


Summary 


e DSP = processing signals with computer algorithms. 
e SSP = statistical DSP = processing in the presence of uncertainties and 
unknowns. 


Review of Linear Algebra 


Vector spaces are the principal object of study in linear algebra. A vector 
space is always defined with respect to a field of scalars. 


Fields 


A field is a set F’ equipped with two operations, addition and mulitplication, 
and containing two special members 0 and 1 (0 ~ 1), such that for all 
{a,b,cleF 


1 l(ea+beF 
2.a+b=b+a 
3.(a+b)+c=a+(b+c) 
4.a+0=a 
5. there exists —a such that a ++ —a = 0 


labe F 

2.ab = ba 

3. (ab)c = a (bc) 

AG. 6 

5. there exists a! such that aa~! = 1 


3.a(b+c)=ab+ac 
More concisely 
1. F is an abelian group under addition 


2. Fis an abelian group under multiplication 
3. multiplication distributes over addition 


Examples 


Q,R,C 


Vector Spaces 


Let F' be a field, and V a set. We say V is a vector space over F if there 
exist two operations, defined for alla € F,uw € V andv € V: 


e vector addition: (u,v) - (w+v) EV 
e scalar multiplication: (a,v) > av € V 


and if there exists an element denoted 0 € V, such that the following hold 
foralac F,b€ Fyanduwe V,veEV,andweVv 


1 lut+(v+w)=(ut+v)+w 
2.uU+v=vU+U 


3.u+0=u 
4. there exists —w such that w+ —u — 0 


More concisely, 


1. V is an abelian group under plus 
2. Natural properties of scalar multiplication 


Examples 


¢ RY” isa vector space over R 

¢ C% is a vector space over C 
Cer R 

° is a vector space over 

¢ RY” is not a vector space over C 


The elements of V are called vectors. 


Euclidean Space 


Throughout this course we will think of a signal as a vector 


The samples {;} could be samples from a finite duration, continuous time 
signal, for example. 


A signal will belong to one of two vector spaces: 


Real Euclidean space 


az € RY (over R) 


Complex Euclidean space 


x € CNX (over C) 


Subspaces 
Let V be a vector space over F’. 


A subset S C V is called a subspace of V if S is a vector space over F' in 
its own right. 


Example: 
V = R?, F=R, S = any line though the origin. 


Ay 


—> X 


A 


S is any line through the origin. 


Are there other subspaces? 


S C V is a subspace if and only if for all a € F and b € F and for all 
s€Sandt€S,(as+bt)eS 


Linear Independence 
Let u1,...,u, € V. 


We say that these vectors are linearly dependent if there exist scalars 
a1,...,@a,% € F such that 


Equation: 
k 
> aQa{,Uy = 0 
i=l 


and at least one a; # 0. 


If [link] only holds for the case ay = ... = az = O,7 we say that the vectors 
are linearly independent. 


Example: 
1 —2 =o 
UP ee Se i) 
2 0 —2 


so these vectors are linearly dependent in R°. 


Spanning Sets 


Consider the subset S = {v1,v2,...,v%}. Define the span of S 


wer 


k 
< S$ >=span(S) = {Soa 
i=l 


Fact: < S > is a subspace of V. 


Example: 


V=R? F-R, S= 101,00},07 = 10 |,0 = |) 1)/> 


< S$ >= xy-plane. 


< S$ > is the xy-plane. 


Aside 


If S is infinite, the notions of linear independence and span are easily 
generalized: 


We say S is linearly independent if, for every finite collection 
U1,---,UK € S, (k arbitrary) we have 


k 
i=1 


The span of Sis 


k 
a 15 = {oa 


dl 


acFAwES A ie<c)| 


Note:In both definitions, we only consider finite sums. 


Bases 


A set B C V is called a basis for V over F if and only if 


1. B is linearly independent 
2-<—b = 


Bases are of fundamental importance in signal processing. They allow us to 


decompose a signal into building blocks (basis vectors) that are often more 
easily understood. 


Example: 
V = (real or complex) Euclidean space, R® or CX. 


B= {ej,...,en} = standard basis 


where the 1 is in the i*" position. 


Example: 
V =CN% over C. 


which is the DFT basis. 


Uk = 


where 2 = VY —l. 


Key Fact 


If Bis a basis for V, then every v € V can be written uniquely (up to order 
of terms) in the form 


where a; € F' and v; € B. 


Other Facts 


e If Sis a linearly independent set, then S' can be extended to a basis. 
e If < S >= V, then S contains a basis. 


Dimension 


Let V be a vector space with basis B. The dimension of V, denoted 
dim (V), is the cardinality of B. 


Every vector space has a basis. 
Every basis for a vector space has the same cardinality. 
= dim (V) is well-defined. 


If dim (V) < oo, we say V is finite dimensional. 


Examples 
vector space field of scalars dimension 
RY R 
Cy C 
ce R 


Every subspace is a vector space, and therefore has its own dimension. 


Example: 
Suppose (S = {u1,...,ux}) C V isa linearly independent set. Then 


came S25) 


Facts 


e If S is a subspace of V, then dim (S) < dim(V). 
e Ifdim(S) = dim(V) < o, thn S=V. 


Direct Sums 
Let V be a vector space, and let S C V and T C V be subspaces. 


We say V is the direct sum of S and 7’, written V = S @ T, if and only if 
for every v © V, there exist unique s € S andt € TJ’ such that v = s+ f. 


If V = S @T, then T is called a complement of S. 


Example: 
V=C' ={f:R—RJfis continuous} 
S = even funcitons inC’ 


T = odd funcitons inC’ 
il i 


BA ee) na ae) Glaciog IN) lilt) 


heap Sevohes eso) Aa oe Seley Ge ee We sie ee ion 
g—g' =h' — his odd and even, which implies g = g’ andh = h’. 


Facts 


1. Every subspace has a complement 
2.V = S@T if and only if 


LSor= {ot 
2.4 5.1 S=V 


3.1fV = S@T, and dim (V) < ov, then 
dim (V) = dim (S$) + dim (T) 


Proofs 


Invoke a basis. 


Norms 


Let V be a vector space over F’. A norm is amapping V — F’, denoted by 
|| + ||, such that forallu €V,v Ee V,andAcCF 


1. || uw ||> Oifu ~0 
2. || Aw ||= |A| || e || 
3. [utes ull+| oll 


Examples 


Euclidean norms: 


2 eR: 


xecn: 


Induced Metric 
Every norm induces a metric on V 
d(u,v) =||u—v | 


which leads to a notion of "distance" between vectors. 


Inner products 


Let V be a vector space over F’, F = Ror C. An inner product is a 
mapping V x V -—> F, denoted (-, -), such that 


1. (v,v) > 0, and (v,v) =O Sv=0 
2. (U,V) = (v,u) 
3. (au + bv, w) = a ((u,w)) + b((v, w)) 


Examples 
RY” over R: 
N 
(x,y) c= (a7 y) = So iy: 
i=l 
CX over C: 


8 
I 


is called the "Hermitian," or "conjugate transpose" of a. 


Triangle Inequality 
If we define || w ||= (w, w), then 
wt v |[<l] ul] + |e | 


Hence, every inner product induces a norm. 


Cauchy-Schwarz Inequality 
ForallweV,ve V, 
(u,v)| <|[ w |] || v || 


In inner product spaces, we have a notion of the angle between two vectors: 


(tu v) = arceos( 7) |) € [0, 2m) 


| w | |e || 


Orthogonality 
uw and v are orthogonal if 

(u,v) =0 
Notation: wu L v. 


If in addition || w |/=|| v ||= 1, we say wu and v are orthonormal. 


In an orthogonal (orthonormal) set, each pair of vectors is orthogonal 
(orthonormal). 


Ay 


= X 


Orthogonal vectors in R?. 


Orthonormal Bases 
An Orthonormal basis is a basis {v; } such that 


lifi=j 


a 
Min Vi) = 949 ree 


Example: 
The standard basis for R% or C” 


Example: 
The normalized DFT basis 


Uk = 


ue 
VN 


Expansion Coefficients 


If the representation of v with respect to {v;} is 
v= » AYU; 
i 
then 
GH Uw) 


Gram-Schmidt 


Every inner product space has an orthonormal basis. Any (countable) basis 
can be made orthogonal by the Gram-Schmidt orthogonalization process. 


Orthogonal Compliments 


Let S C V bea subspace. The orthogonal compliment S is 


St={ulueV A ((u,v) =0) A Vu: (ve S)} 
S'+ is easily seen to be a subspace. 


If dim (v) < 00, then V= SQ S". 


Note:If dim (v) = 00, then in order to have V = S$ @ S+ we require V to 
be a Hilbert Space. 


Linear Transformations 


Loosely speaking, a linear transformation is a mapping from one vector 
space to another that preserves vector space operations. 


More precisely, let V, W be vector spaces over the same field F’. A linear 
transformation is a mapping 7’: V — W such that 


T(au + bv) = aT(u) + bT(v) 
foralac F,b€ Fanduc V,veEV. 


In this class we will be concerned with linear transformations between (real 
or complex) Euclidean spaces, or subspaces thereof. 


Image 


image (T) = {w| we W A T(v) = wfor somev} 


Nullspace 


Also known as the kernel: 


ker (T) = {vj v EV A (T(v) = 90)} 


Both the image and the nullspace are easily seen to be subspaces. 


Rank 


rank (T) = dim (image (T)) 


Nullity 


null (7) = dim (ker (T')) 


Rank plus nullity theorem 


rank (TJ) + null (7) = dim (V) 


Matrices 


Every linear transformation T' has a matrix representation. If 
T:EX ~ E”“,E=RorC, then T is represented by an M x N matrix 


Q@M1 --- QA@MN 


where (aj;...ayi)’ = T(e;) and e; = (0...1...0)" is the i** standard 
basis vector. 


Note:A linear transformation can be represented with respect to any bases 
of EN and E™, leading to a different A. We will always represent a linear 
transformation using the standard bases. 


Column span 


colspan (A) =< A >=image (A) 


Duality 


If A: RY > R™, then 


ker+(A) =image (A*) 


HA<:<C* —-C™ , then 


ker~ (A) =image (A®) 


Inverses 


The linear transformation/matrix A is invertible if and only if there exists a 
matrix B such that AB = BA = I (identity). 


Only square matrices can be invertible. 
Let A: FX — F% be linear, F = R or C. The following are equivalent: 
1. A is invertible (nonsingular) 


2.rank(A) = N 
3. null(A) = 0 


4.det A #0 
5. The columns of A form a basis. 


If A~' = A? (or A# in the complex case), we say A is orthogonal (or 
unitary). 


The Bayesian Paradigm 


Statistical analysis is fundamentally an inversion process. The objective is to the 
"causes'"--parameters of the probabilistic data generation model--from the 
"effects"--observations. This can be seen in our interpretation of the likelihood 
function. 


Given a parameter @, observations are generated according to 
p (| @) 


The likelihood function has the same form as the conditional density function 
above 


1(O|x) =p (x/@) 


except now 2 is given (we take measurements) and @ is the variable. The 
likelihood function essentially inverts the role of observation (effect) and 
parameter (cause). 


Unfortunately, the likelihood function does not provide a formal framework for the 
desired inversion. 


One problem is that the parameter 0 is supposed to be a fixed and deterministic 
quantity while the observation z is the realization of a random process. So their 
role aren't really interchangeable in this setting. 


Moreover, while it is tempting to interpret the likelihood /(@|x) as a density 
function for 9, this is not always possible; for example, often 


[ ee) ae 


Another problematic issue is the mathematical formalization of statements like: 
"Based on the measurements 2, I am 95% confident that @ falls in a certain range." 


Example: 

Suppose you toss a coin 10 times and each time it comes up "heads." It might be 
reasonable to say that we are 99% sure that the coin is unfair, biased towards 
heads. 


Formally: 
Ho : 0 = prob heads > 0.5 


ze en gu*(1 9)" >* 


which is the binomial likelihood. 
pte > 0.5 |e) —<7 
The problem with this is that 
p (6 € Holz) 


implies that 0 is a random, not deterministic, quantity. So, while "confidence" 
statements are very reasonable and in fact a normal part of "everyday thinking," 
this idea can not be supported from the classical perspective. 


All of these "deficiencies" can be circumvented by a change in how we view the 
parameter 0. 


If we view @ as the realization of a random variable with density p (0), then Bayes 
Rule (Bayes, 1763) shows that 


— p(w/6)p (6) 
BOI?) =" (a|8) p (0) d 8 


Thus, from this perspective we obtain a well-defined inversion: Given 2, the 
parameter @ is distributing according to p (0| a). 


From here, confidence measures such as p (9 € Ho|a) are perfectly legitimate 
quantities to ask for. 


Bayesian statistical model 
A statistical model compose of a data generation model, p (a|@), and a prior 
distribution on the parameters, p (8). 


The prior distriubtion (or prior for short) models the uncertainty in the 
parameter. More specifically, p (8) models our knowledge--or lack thereof--prior 


to collecting data. 
Notice that 


p (#|) p (8) 


p (O|z) = a 


xp (|) p (8) 


since the data a are known, p (2) is just a constant. Hence, p (0| x) is 
proportional to the likelihood function multiplied by the prior. 


Bayesian analysis has some significant advantages over classical statistical 
analysis: 


1. properly inverts the relationship between causes and effects 

2. permits meaningful assessments in confidence regions 

3. enables the incorporation of prior knowledge into the analysis (which could 
come from previous experiments, for example) 

4. leads to more accurate estimators (provided the prior knowledge is accurate) 

5. obeys the Likelihood and Sufficiency principles 


Example: 
Vigta=41 Ns (te Ae WV) 
Wy, ~ V0, 02) 
iid. 
Be atlnne 
Ap W dn 


MVUB and MLE estimator. Now suppose that we have prior knowledge that 
—Ag < A < Apo. We might incorporate this by forming a new estimator 
Equation: 
—Apo if A < —Ap 
A= A if aA wt A 
Ao if A > Ao 


This is called a truncated sample mean estimator of A. Is A a better estimator of 
A than the sample mean A? 

Let p (a) denote the density of A. Since A = = Gy oe = NA, ¢ J. 
The density of Ais given by 

Equation: 


Sg) = Pe [A z —Ao| SG AE (Olin IDs [A > Ay| 5(a — Ag) 


p{a) 


a 


. -Ap Ag 
Now consider the MSE of the sample mean A. 
Equation: 
De, [e@) 
MSE (A) =| (a— A)’ p(a)da 
—Co 
Note 


1. A is biased ([link]). 
2. Although A is MVUB, A is better in the MSE sense. 
3. Prior information is aptly described by regarding A as a random variable with 


a prior distribution U(— Ap, Ao), which implies that we know 
—Apy < A < Ap, but otherwise A is abitrary. 


Mean of A = A. Mean of A # A. 


p(a) 


-Ag A a * 


fab) 


randomly 
generate 
8 


pe eats cia ee re rere eal 


Where w is the noise and x is the observation. 


Example: 


Wie ye lee SN peat eA 


Prior distribution allows us to incorporate prior information regarding unknown 
paremter--probable values of parameter are supported by prior. Basically, the prior 


reflects what we believe "Nature" will probably throw at us. 


Elements of Bayesian Analysis 


¢ (a) joint distribution 
p (x, 8) =p (x|8) p (8) 


e (b) marginal distributions 


where p (8) is a prior. 
¢ (c) posterior distribution 


4 i 
p (6|#) = > = 
ple 


Example: 


v8,8-¢ [0,1]: (p (w|9) = ("ora —4)"~*] 


which is the Binomial likelihood. 


1 


= p-1 
p (6) = = 6"1(1 — 6) 
B(a, B) 
which is the Beta prior distriubtion and B(a, 8) = Gee 
p(9) 
0 a 1 6 
at+ 8 


This reflects prior 
knowledge that 
most probable 
values of a are 


close to aaa 


Joint Density 


ne 
L 

sa NE pane Hy n—2£+B-1 
B(a, B) - 


marginal density 


= (") (a+) I(a+a)P(n—-x+8) 
P T(ar(8) I(atB+n) 


posterior density 


Qota—1lgB+n—a-1 


os a aerate x) 


where B(a + x, 8 +n — 2) is the Beta density with parameters 
a’ =a+azandf’=B+n-2z 


Selecting an Informative Prior 


Clearly, the most important objective is to choose the prior p (0) that best reflects 
the prior knowledge available to us. In general, however, our prior knowledge is 
imprecise and any number of prior densities may aptly capture this information. 
Moreover, usually the optimal estimator can't be obtained in closed-form. 


Therefore, sometimes it is desirable to choose a prior density that models prior 
knowledge and is nicely matched in functional form to p (a|@) so that the optimal 
esitmator (and posterior density) can be expressed in a simple fashion. 


Choosing a Prior 
1. Informative Priors 


e design/choose priors that are compatible with prior knowledge of unknown 
parameters 


2. Non-informative Priors 


e attempt to remove subjectiveness from Bayesian procedures 
e designs are often based on invariance arguments 


Example: 
Suppose we want to estimate the variance of a process, incorporating a prior that 
is amplitude-scale invariant (so that we are invariant to arbitrary amplitude 


rescaling of data). p (s) = + satisifies this condition. 0? ~p (s) = Ao? ~p (s) 
where p (s) is non-informative since it is invariant to amplitude-scale. 


Conjugate Priors 


Idea 


Given p (a|6), choose p (0) so that p (6|a) xp (a|0) p (0) has a simple 
functional form. 


Conjugate Priors 


Choose p (0) € #, where P is a family of densities (e.g., Gaussian family) so 
that the posterior density also belongs to that family. 


conjugate prior 
p (6) is a conjugate prior for p (a|0) ifp (0) € A=>p (Ola) ce FP 


Example: 
Vinge Ala wn aaiN (nee — Ay) 
W,~ AN (0,5) 


iid. Rather than modeling A ~ U(—Ap, Ao) (which did not yield a closed-form 
estimator) consider 


p(A) 


} A 


With w = 0 anda, = + Ao this Gaussian prior also reflects prior knowledge that 
it is unlikely for |A| > Apo. 


The Gaussian prior is also conjugate to the Gaussian likelihood 


p (@| A) = bet Bhar a 


so that the resulting posterior density is also a simple Gaussian, as shown next. 


First note that 


= N =1 2 7 
p (x|A) = = ee tn ae (NA 2NAz) 
(2707) 2 
ms N 
where x = 4 pan ae 
Equation: 
A)p(A) 
A a“ = __ p(z|A)p(A) _ 
pe) Jp(@| A)p(A)dA 
=1/ 1 (yy2_onaz ee ?) 
(a cals )+ shat ) 
S1( 4 .(wa2-2Nn az )+—+,(4-p)? 
f2e 2 (2 (wa 2NA )+ aru 1) ed 
= Sa (A) 
f~ era 4 
where Q(A) = = A ae a aan ca Now let 


il 
CAle = 
2 = be 
ea (3 , iz) ae 
Then by "completing the square" we have 
Equation: 
Hala” : 
Q(A) = gir (A? -2najeA + wale”) — Fe + $a 
= i 2 HAle be 
7 CAlx? (A 7 HAle) - TAlx* i oA 
Hence, 


where Math input error is the "unnormalized" Gaussian density and 


2 2 
= & — fae, | is a constant, independent of A. This implies that 
1 re (A-pale)” 
p (Ala) = eA 
270 Ale 2 


where Ala ~ N (Males The"). Now 
Equation: 
A = E[Ala] 
= {[Ap(Ala)dA 
= KHAle 


Where 0 < a = aE x1 


oaet a 
Interpretation 
1. When there is little datag 4? < x qa is small and A oe 
2. When there is a lot of data 042 >> an a~landA=z. 


Interplay Between Data and Prior Knowledge 
Small N — A favors prior. 


Large N — A favors data. 


The Multivariate Gaussian Model 


The multivariate Gaussian model is the most important Bayesian tool in signal 
processing. It leads directly to the celebrated Wiener and Kalman filters. 


Assume that we are dealing with random vectors z and y. We will regard y as a 
signal vector that is to be estimated from an observation vector 2. 


y plays the same role as @ did in earlier discussions. We will assume that y is px1 
and x is Nx1. Furthermore, assume that x and y are jointly Gaussian distributed 


(AG ea) 


Ela] = 0, Ely] = 0, Elwa"] = Ryx, Elwy"| = Rey, Elyx"| = Ryx, 
Ri 
Re XX > 
(z Ryy 


Example: 
2=y+Ww,W~Y(0,07l) 


P (y) = V(0, Ryy) 


which is independent of W. E|a| = Ely] + E|W] =0, 
Elxx"| = Elyy’ | + E|yw*| + E|Wy’| + E|Ww"*|] = Ryy + ot, 
Eley") = Elyy’| + E|Wy"] = Ry = Bye’). 


)~(() a, ae) 
y 0 Ryy Ryy 
From our Bayesian perpsective, we are interested in p (y|a). 
Equation: 


p(x,y) 


p(x) 


p (y|a) 
(2m) (2n)-3 (det Ryte? eae (5) 


N —-1 —-1 = 
(2n)~ 2 (det Rex) Ze 2” B® 


In this formula we are faced with 


ae ts sel ae 
Ryx Ryy 


The inverse of this covariance matrix can be written as 


feels eet 


where Q = Ry, — RyxRyxRyy. (Verify this formula by applying the right hand 
side above to R to get I.) 


Sufficient Statistics 


Introduction 


Sufficient statistics arise in nearly every aspect of statistical inference. It is 
important to understand them before progressing to areas such as 
hypothesis testing and parameter estimation. 


Suppose we observe an V-dimensional random vector X, characterized by 
the density or mass function fg (a), where @ is a p-dimensional vector of 
parameters to be estimated. The functional form of f (a) is assumed known. 
The parameter 90 completely determines the distribution of X. Conversely, 
a measurement a of X provides information about @ through the 
probability law fg (2). 


Example: 


xX 
Suppose X = ey , where X; ~ V(0, 1) are IID. Here @ is a scalar 
2 


parameter specifying the mean. The distribution of X is determined by 0 
through the density 


1 ce eee __ (22-0)? 
fy (x) = ——e 
/ 24 V/ 20 


100 
On the other hand, if we observe x = ( 


rel , then we may safely assume 


0 = 0 is highly unlikely. 


The N-dimensional observation X carries information about the p- 
dimensional parameter vector @. If p < N, one may ask the following 
question: Can we compress 2 into a low-dimensional statistic without any 
loss of information? Does there exist some function £ = Ta), where the 


dimension of t is MZ < NN, such that t carries all the useful information 
about 0? 


If so, for the purpose of studying @ we could discard the raw measurements 
x and retain only the low-dimensional statistic £. We call ¢ a sufficient 
statistic. The following definition captures this notion precisely: 


Let Xj,..., X.¢ be arandom sample, governed by the density or 

probability mass function f (a |@). The statistic Ta) is sufficient for 
6 if the conditional distribution of x, given T(a) = t, is independent 
of @. Equivalently, the functional form of fg), (a) does not involve 0. 


How should we interpret this definition? Here are some possibilities: 


1. Let fg (a, €) denote the joint density or probability mass function on 
(X,7T(X)). If T(X) is a sufficient statistic for 8, then 
Equation: 


fo (a) = fo (x, T(x)) 


fo\z (x) fo (t) 
= f(a|t) fo (t) 


Therefore, the parametrization of the probability law for the measurement z 
is manifested in the parametrization of the probability law for the statistic 
Ta), 


2. Given t = T(a), full knowledge of the measurement 2 brings no 
additional information about @. Thus, we may discard a and retain on the 
compressed statistic f. 


3. Any inference strategy based on fg (a) may be replaced by a strategy 
based on fg (t). 


Example: 
Binary Information Source 


(Scharf, pp.78) Suppose a binary information source emits a 
sequence of binary (0 or 1) valued, independent variables 
£1,---,xN. Each binary symbol may be viewed as a 
realization of a Bernoulli trial: z,, ~ Bernoulli (6), iid. The 
parameter 6 € [0, 1] is to be estimated. 


The probability mass function for the random sample 
bs (ay. a ay)" is 


Equation: 
N N 
fg (x) = || fo (en) [[ *G-9)*" 
n=1 n=1 


where k = ee Xp, is the number of 1's in the sample. 
We will show that k is a sufficient statistic for x. This will 
entail showing that the conditional probability mass 
function fg), (a) does not depend on 8. 


The distribution of the number of ones in NV independent 
Bernoulli trials is binomial: 


fy (k) = (.) oy = 


Next, consider the joint distribution of (a, 5) 2,). We have 


Hy ()) (« Ys] 


Thus, the conditional probability may be written 
Equation: 


fol, (@) = ne 


This shows that k is indeed a sufficient statistic for 0. The 
N values 21,...,y can be replaced by the quantity k 
without losing information about 0. 


Exercise: 
Problem: 
In the previous example, suppose we wish to store in memory the 


information we possess about 8. Compare the savings, in terms of bits, 
we gain by storing the sufficient statistic k instead of the full sample 


L1,+++ UN. 
Determining Sufficient Statistics 


In the example above, we had to guess the sufficient statistic, and work out 
the conditional probability by hand. In general, this will be a tedious way to 


go about finding sufficient statistics. Fortunately, spotting sufficient 
statistics can be made easier by the Fisher-Neyman Factorization Theorem. 


Uses of Sufficient Statistics 


Sufficient statistics have many uses in statistical inference problems. In 
hypothesis testing, the Likelihood Ratio Test can often be reduced to a 
sufficient statistic of the data. In parameter estimation, the Minimum 
Variance Unbiased Estimator of a parameter @ can be characterized by 
sufficient statistics and the Rao-Blackwell Theorem. 


Minimality and Completeness 


Minimal sufficient statistics are, roughly speaking, sufficient statistics that 
cannot be compressed any more without losing information about the 
unknown parameter. Completeness is a technical characterization of 
sufficient statistics that allows one to prove minimality. These topics are 
covered in detail in this module. 


Further examples of sufficient statistics may be found in the module on the 
Fisher-Neyman Factorization Theorem. 


The Fisher-Neyman Factorization Theorem 


Determining a sufficient statistic directly from the definition can be a 
tedious process. The following result can simplify this process by allowing 
one to spot a sufficient statistic directly from the functional form of the 
density or mass function. 

Fisher-Neyman Factorization Theorem 


Let fg (a) be the density or mass function for the random vector a, 
parametrized by the vector @. The statistic ¢ = Ta) is sufficient for 0 if 
and only if there exist functions a(a) (not depending on 0) and bg (€) such 
that 


fo (a) = a(@) bo (¢) 


for all possible values of x. 

In an earlier example we computed a sufficient statistic for a binary 
communication source (independent Bernoulli trials) from the definition. 
Using the above result, this task becomes substantially easier. 


Example: 
Bernoulli Trials Revisited 
Suppose z,, ~ Bernoulli (8) are IID, 
Vit — peer dee) Denote 


cel ae De Then 
Equation: 


fy (@) = Tp O?(1- 0) 
K=O). a 
a(a) bg (k) 


where k = Soa foe) — land 
bg (k) = 6*(1 — 6)” “. By the Fisher-Neyman 
factorization theorem, k is sufficient for 0. 


The next example illustrates the appliction of the theorem to a continuous 
random variable. 


Example: 
Normal Data with Unknown Mean 


Consider a normally distributed random sample 
£1,---)£N ~ WV(6,1), IID, where @ is unknown. The joint 
pdf of # = (21...ap)" is 


We would like to rewrite fg (a) is the form of a(a) bg (€), 
where dim (t) < N. At this point we require a trick-one 
that is commonly used when manipulating normal densities, 
and worth remembering. Define z = + Sey Ln, the 
sample mean. Then 

Equation: 


fg (a) = (+) et 


aN 
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Now observe 
Equation: 


so the middle term vanishes. We are left with 


N 


fa (#) = (=) ? PEL (2) FEL) 


- NO?) 
where a(x) = (3~) F et Lint (2-2) 
3 rhe (2-4) 
Bat hier. , and t = a. Thus, the sample 


mean is a one-dimensional sufficient statistic for the mean. 


Proof of Theorem 


First, suppose £ = T(x) is sufficient for @. By definition, fg) 7(~)-¢ (a) is 
independent of @. Let fg (a, t) denote the joint density or mass function for 
(X, T'(X)). Observe fg (a) =fo (a, ¢). Then 

Equation: 


fo (a) = fo (x,t) 
fo\z (x) fo (t) 
a(a) be (t) 


where a(a) =fg)¢ (a) and be (t) =f (£). We prove the reverse 
implication for the discrete case only. The continuous case follows a similar 


pp.127). 
Suppose the probability mass function for a can be written 
fg (w) = ala) bo (x) 


where t = Ta). The probability mass function for ¢ is obtained by 
summing fg (a, ¢) over all a such that T(a) = t: 
Equation: 


fg (t) = DUT (a)=t fo 


| 
p= 
S 
I 
~~ 
3” 


Therefore, the conditional mass function of z, given f, is 
Equation: 


fox (a) = 08) 


yr(e)=t a(a) 


This last expression does not depend on @, so t is a sufficient statistic for 0. 
This completes the proof. 


Note:From the proof, the Fisher-Neyman factorization gives us a formula 
for the conditional probability of a given €. In the discrete case we have 


a(x) 


DT (@)=t a(x) 


An analogous formula holds for continuous random variables (Scharf, 
pp-82). 


et 


Further Examples 


The following exercises provide additional examples where the Fisher- 
Neyman factorization may be used to identify sufficient statistics. 
Exercise: 


Problem: 
Uniform Measurements 


Suppose 21,...,y are independent and uniformly distributed on the 
interval [0,, 05]. Find a sufficient statistic for 9 = (0,02)*. 


Note:Express the likelihood fg (a) in terms of indicator functions. 


Exercise: 
Problem: 
Poisson 


Suppose 21,..., 2 are independent measurements of a Poisson 
random variable with intensity parameter 0: 


—092 
Wie = 051 2? (1 (2) = mt ) 


x! 


Find a sufficient statistic ¢ for 0. 
What is the conditional probability mass function of a, given t, where 
P= (Pi. ay)? 
Exercise: 
Problem: 


Normal with Unknown Mean and Variance 


Consider x1,...,2N ~ Ai, 0”), IID, where 6; = ps and 62 = 0? 
are both unknown. Find a sufficient statistic for 9 = (0,02)°. 


Note: Use the same trick as in [link]. 


Signal Classifications and Properties 
Describes various classifications of signals. 


Introduction 


This module will begin our study of signals and systems by laying out some 
of the fundamentals of signal classification. It is essentially an introduction 
to the important definitions and properties that are fundamental to the 
discussion of signals and systems, with a brief discussion of each. 


Classifications of Signals 


Continuous-Time vs. Discrete-Time 


As the names suggest, this classification is determined by whether or not 
the time axis is discrete (countable) or continuous ((link]). A continuous- 
time signal will contain a value for all real numbers along the time axis. In 
contrast to this, a discrete-time signal, often created by sampling a 
continuous signal, will only have values at equally spaced intervals along 
the time axis. 


this axis continuous 
or discrete 


Analog vs. Digital 


The difference between analog and digital is similar to the difference 
between continuous-time and discrete-time. However, in this case the 
difference involves the values of the function. Analog corresponds to a 


continuous set of possible function values, while digital corresponds to a 
discrete set of possible function values. An common example of a digital 
signal is a binary sequence, where the values of the function can only be 
one or zero. 


this axis continuous 
or discrete 


Periodic vs. Aperiodic 


Periodic signals repeat with some period 7’, while aperiodic, or 
nonperiodic, signals do not ([{link]). We can define a periodic function 
through the following mathematical expression, where ¢ can be any number 
and T is a positive constant: 

Equation: 


f(t) = f(t +T) 


fundamental period of our function, f(t), is the smallest value of T that 
the still allows [link] to be true. 


A periodic signal with period To 


An aperiodic signal 


Finite vs. Infinite Length 


Another way of classifying a signal is in terms of its length along its time 
axis. Is the signal defined for all possible values of time, or for only certain 
values of time? Mathematically speaking, f(t) is a finite-length signal if it 
is defined only over a finite interval 


= t <1 


where t; < tg. Similarly, an infinite-length signal, f(t), is defined for all 
values: 


—o <t< oo 


Causal vs. Anticausal vs. Noncausal 


Causal signals are signals that are zero for all negative time, while 
anticausal are signals that are zero for all positive time. Noncausal signals 
are signals that have nonzero values in both positive and negative time 


(Link). 


f(t) 


zero here 


A causal signal 


f(t) 


zero here 


An anticausal signal 


f(t) 


A noncausal signal 


Even vs. Odd 


An even signal is any signal f such that f(t) = f(—t). Even signals can be 
easily spotted as they are symmetric around the vertical axis. An odd 
signal, on the other hand, is a signal f such that f(t) = — f(—t) ([link)). 


fe(t) 
es. t 


An even signal 


fo(t) 


An odd signal 


Using the definitions of even and odd signals, we can show that any signal 
can be written as a combination of an even and odd signal. That is, every 
signal has an odd-even decomposition. To demonstrate this, we have to look 
no further than a single equation. 

Equation: 


fl) = 5 (F®) + FH) + 5 FO — FH) 


By multiplying and adding this expression out, it can be shown to be true. 
Also, it can be shown that f(t) + f(—t) fulfills the requirement of an even 
function, while f(t) — f(—t) fulfills the requirement of an odd function 
({link]). 


Example: 


The signal we will 
decompose using odd- 
even decomposition 


f(t) 


f(t) : 
“1 


Even part: e(t) = + (f(t) + f(—t)) 


oe 


2e(t) 
2 
t 
4 


it's even! 


- 


f(t) 
2o(t) 


oe 


it's odd! 


oe 


Odd part: o(t) = + (f(t) — f(t) 


Check: e(¢) + o(t) = f(t) 


Deterministic vs. Random 


A deterministic signal is a signal in which each value of the signal is fixed, 
being determined by a mathematical expression, rule, or table. On the other 
hand, the values of a random signal are not strictly defined, but are subject 
to some amount of variability. 


Deterministic Signal 


0 sf, ype yoo way bi TON 2 


Random Signal 


Example: 
Consider the signal defined for all real ¢ described by 
Equation: 


_ fsin (2rt)/t t>1 
ro={ 0 peel 


This signal is continuous time, analog, aperiodic, infinite length, causal, 
neither even nor odd, and, by definition, deterministic. 


Signal Classifications Summary 


This module describes just some of the many ways in which signals can be 
classified. They can be continuous time or discrete time, analog or digital, 
periodic or aperiodic, finite or infinite, and deterministic or random. We can 
also divide them based on their causality and symmetry properties. 


Linear Models 


Finding an MVUB estimator is a very difficult task, in general. However, a 
large number of signal processing problems can be respresented by a linear 
model of the data. 


Importance of Class of Linear Models 


1. MVUB estimator within this class is immediately evident 
2. Statistical performance analysis of linear models is very 
straightforward 


General Form of Linear Model (LM) 
«= HO+w 


where a is the observation vector, H is the known matrix (observation or 
system matrix), 8 is the unknown parameter vector, and w is the vector 
of White Guassian noise w ~ (0,07). 


Example: 


Vn metlecs N} a(t, = ASB, ew) 


Ly 
4 = 

LN 

Wi 
Ww — 


WN 


H= 
1 N 
Probability Model for LM 
a= HO+w 


a ~p (z|0) = (HO, 071) 
CRLB and NVUB Estimator 


@ = g(a) the MVUB estimator iff 


Olog p (2| 9) 


Pacis = 1(8) (9(«) - 8) 


In the case of the LM, 


dlog p (x| 0) 1 \ 0 (a? a — 2x7 HO + 67H" HO) 
a0 ( ) a0 


Now using identities 


8 (b76) 
es 
Q(0TAO) ve 


for A symmetric. 
We have 
O01 0 1 
OREO ENO (cH) 
06 o? 


Assuming H? H is invertible 


00 o? 
which leads to 
Equation: 
MVUB Estimator 
6=(H7H) HTe 
Equation: 


Fisher Information Matrix 


_ HH 


o2 


1() 


MVUB Estimator for the LM 
If the observed data can be modeled as 
z= HO+w 


where w ~ VY (0, oI ) and # is invertible. Then, the MVUB estimator is 


and the covariance of @ is 


and @ attains the CRLB. 
Note: 6 ~ N(8, o”(H" H) 2) 


Linear Model Examples 


Example: 

Curve Fitting 
[missing resource: | 
Model: 


Vn,n € {1,...,.N}: (a(tp) = 01 + Ooty +... + Opty? | + w(tn)) 


where 0; + O9tn +... + Oytn?} is a (p — 1)**-order polynomial and 
w(tn) ~ V (0, o”) idd. Therefore, 


z— HO+w 
x(t) 


A, 
cal (ee 
A» 
1S ai foe 
1 to ee 
H= 
Stat ee ty?! 


where H is the Vandermonde matrix. The MVUB estimator for @ is 


il 


ee ee saat: 
Example: 
System Identification 
[missing_resource: | 
m—1 
He) hlk]z* 
k=0 


Vn, n € {0,...,N—1}: (on = 


Where w[n] ~ 4 (0,07) idd. Given « and u, estimate h. 
In matrix form 


u(1] u(0] 0 
t= + Ww 
ulN—1] w[N—2) ... ulN—m)) \nIw—m 

where 

u|0) 0 0 

uf! u(0] 0 

; = sp 
u[N — 1] Ane 2] u[N — ml] 
and 
h{0| 
3 = 0 
h|.N Z m| 
Equation: 
MVUB estimator 


6=(H"H) ‘H' 


(3°) =o (8) =, 


An important question in system identification is how to choose the input 
un] to "probe" the system most efficiently. 
First note that 


o(4;) ° = ei; Cze; 


where e; = (0...010...0)”. Also, since C’ 57 is symmetric positive 
definite, we can factor it by 


Co} 


where D is invertible.[ footnote] Note that 
Equation: 


(7D? (DT) *e,) =1 


The Schwarz inequality shows that [link] can become 
Equation: 


1 < (e;"D* De;) (ga (D*) “Ss) 


ce (e:7C; *e:) (e:* C;e:) 

which leads to 

Pee il 2 

eC,” ej (H A); ; 
The minimum variance is achieved when equality is attained in [link]. This 
happens only if 7; = De; is proportional to 72 = D? e;. That is, 
71 = C'n2 for some constant C’. Equivalently, 
Mid © le een (D? De; = Cen) 

Ho 


o2 


iD) = Ca — 


which leads to 


fees es 


o2 


es Ee 


Combining these equations in matrix form 


Ci Oe sze 0 
Ht'H=07|0 ce... O 
Ope Tae eene ss 


Therefore, in order to minimize the variance of the MVUB estimator, u|n| 
should be chosen to make H? H diagonal. 
Cholesky factorization 


ae jie 5 ag ol eae 100) (ati), = Yuin ~sae 


For large N, this can be approximated to 
N13 
(H7H),.~ SY > ulnjufn+li—3il 


tj 
m—0 


using the autocorrelation of seq. u|n]. 


Note: u[n] = 0 forn < Oandn > N — 1, letting limit of sum —co, oo 
gives approx. 


These steps lead to 


Taal Taalv 
Ny [1] [0] 
Tuulm—1) ruulm — 2] Ran 
Pell Tall Peli — 1) 
eal ere 01 a : 
where is the Toeplitz 


Temnl0l 


— 


autocorrelation matrix and 


For H* H to be diagonal, we require r[k] = 0, k 4 0. This condition is 
approximately realized if we take u[n| to be a pseudorandom noise 
sequence (PRN)[footnote]. Furthermore, the PRN sequence simplifies the 
estimator Computation: 


which leads to 


where ae u[n — i]2[n] = Nrux[i]. rux|2] is the cross-correlation 


between input and output sequences. 
maximal length sequences 


Hence, the approximate MVUB estimator for large N with a PRN input is 


Vi,i € {0,1,....m—1}: (‘i = Pall ) 


anol 
il N-1-1 
fee | = FF a u|n|z[n + 4] 
7; we : 
errr a oe u‘|n| 


CRLB for Signal in White Gaussian Noise 


Vit @ {lieu Ns (ty = 8a 8) wa) 


1 (© S* (es (6 
p (e|6) = 2 eax Ea (a 20(0))) 


Olog p (x|8@) 1 s 


Hypothesis Testing 


Suppose you measure a collection of scalars x;,..., 2. You believe the 
data is distributed in one of two ways. Your first model, call it Ho, 
postulates the data to be governed by the density f(x) (some fixed 
density). Your second model, Hj, postulates a different density f1(z). 
These models, termed hypotheses, are denoted as follows: 


Hy itn Jole) n= 1LacN 
Hitt fie) n= le 


A hypothesis test is a rule that, given a measurement z, makes a decision 
as to which hypothesis best "explains" the data. 


Example: 

Suppose you are confident that your data is normally distributed with 
variance 1, but you are uncertain about the sign of the mean. You might 
postulate 


Ag ty, 7 (1, 1) 
ie ey (a) 


These densities are depicted in [link]. 


Assuming each hypothesis is a priori equally likely, an intuitively 
appealing hypothesis test is to compute the sample mean x = + SS a Ba 


, and choose Ho if x < 0, and H; if x > 0. As we will see later, this test is 
in fact optimal under certain assumptions. 


Generalizations and Nomenclature 


The concepts introduced above can be extended in several ways. In what 
follows we provide more rigorous definitions, describe different kinds of 
hypothesis testing, and introduce terminology. 


Data 


In the most general setup, the observation is a collection 21,..., xy of 
random vectors. A common assumption, which facilitates analysis, is that 
the data are independent and identically distributed (IID). The random 
vectors may be continuous, discrete, or in some cases mixed. It is generally 
assumed that all of the data is available at once, although for some 
applications, such as Sequential Hypothesis Testing, the data is a never 
ending stream. 


Binary Versus M-ary Tests 


When there are two competing hypotheses, we refer to a binary hypothesis 
test. When the number of hypotheses is M > 2, we refer to an M-ary 
hypothesis test. Clearly, binary is a special case of /-ary, but binary tests 
are accorded a special status for certain reasons. These include their 
simplicity, their prevalence in applications, and theoretical results that do 
not carry over to the M-ary case. 


Example: 
Phase-Shift Keying 


Suppose we wish to transmit a binary string of length r over 
a noisy communication channel. We assign each of the 

M = 2" possible bit sequences to a signal s*, 

k = {1,...,M} where 


2m(k— I 
s* = cos (2nfon ale ae) 


This symboling scheme is known as phase-shift keying 
(PSK). After transmitting a signal across the noisy channel, 
the receiver faces an M-ary hypothesis testing problem: 


H):x2=s'+w 


Hy :2=s@+w 


where w ~ N (0, ca 


In many binary hypothesis tests, one hypothesis represents the absence of a 
ceratin feature. In such cases, the hypothesis is usually labelled Ho and 
called the null hypothesis. The other hypothesis is labelled Hy and called 
the alternative hypothesis. 


Example: 
Waveform Detection 


Consider the problem of detecting a known signal 
s = (s1...8y)’ in additive white Gaussian noise (AWGN). 
This scenario is common in sonar and radar systems. 
Denoting the data as @ = (x1...ay)", our hypothesis 
testing problem is 
dice ont} 

A, :x2=s+w 
where w ~ WV (0, od Hg is the null hypothesis, 
corresponding to the absence of a signal. 


Tests and Decision Regions 


Consider the general hypothesis testing problem where we have WN d- 
dimensional observations 21,...,£j and M hypotheses. If the data are 
real-valued, for example, then a hypothesis test is a mapping 


y: (R*)” > {1,...,M} 


For every possible realization of the input, the test outputs a hypothesis. 
The test y partitions the input space into a disjoint collection R,,..., Ry, 
where 


Re = A Giges Qn) Ol Cisse stn SF} 


The sets Rx are called decision regions. The boundary between two 
decision regions is a decision boundary. [link] depicts these concepts when 
d=2N = 1, ond = 3. 


Simple Versus Composite Hypotheses 


If the distribution of the data under a certain hypothesis is fully known, we 
call it a simple hypothesis. All of the hypotheses in the examples above are 
simple. In many cases, however, we only know the distribution up to certain 
unknown parameters. For example, in a Gaussian noise model we may not 
know the variance of the noise. In this case, a hypothesis is said to be 
composite. 


Example: 
Consider the problem of detecting the signal 


Si = COS ota ao) Vl lee 
where & is an unknown delay parameter. Then 
el Gye 
A,:x2=s+w 


is a binary test of a simple hypothesis (Hg) versus a composite alternative. 
Here we are assuming w, ~ WV (0, an, with ao? known. 


Often a test involving a composite hypothesis has the form 
A QO: d= Ao 
H,:04% 


where 6p is fixed. Such problems are called two-sided because the 
composite alternative "lies on both sides of Hg." When @ is a scalar, the test 


Ho :9< 4% 
H,:0> 4 


is called one-sided. Here, both hypotheses are composite. 


Example: 

Suppose a coin turns up heads with probability p. We want to assess 
whether the coin is fair (p = 5). We toss the coin NV times and record 
L1,--+-;LN (Lp = 1 means heads and z,, = O means tails). Then 


is a binary test of a simple hypothesis (Hg) versus a composite alternative. 
This is also a two-sided test. 


Errors and Probabilities 


In binary hypothesis testing, assuming at least one of the two models does 
indeed correspond to reality, there are four possible scenarios: 


e Case 1 Hp is true, and we declare Ho to be true 
e Case 2 Hy is true, but we declare H; to be true 
e Case 3 H; is true, and we declare H to be true 
e Case 4 H; is true, but we declare Ho to be true 


In cases 2 and 4, errors occur. The names given to these errors depend on 
the area of application. In statistics, they are called type I and type II 
errors respectively, while in signal processing they are known as a false 
alarm or a miss. 


Consider the general binary hypothesis testing problem 
Ho :@%~ fo(a),0 E Oo 
Ai: a2wrw fo(a), 0 EO; 


If Ho is simple, that is, 09 = {6}, then the size (denoted a), also called 
the false-alarm probability (Pr), is defined to be 


a = Pr = Pr[6o, declareH}| 
When @, is composite, we define 
a = Pr = supgce, (Pr/0, declareH;]) 


For 8 € @,, the power (denoted {), or detection probability (Pp), is 
defined to be 


B = Pp = Pr{6, declareH,| 
The probability of a type II error, also called the miss probability, is 
Py=1-—-Pp 


If H; is composite, then 8 = 6(6) is viewed as a function of 6. 


Criteria in Hypothesis Testing 


The design of a hypothesis test/detector often involves constructing the 
solution to an optimization problem. The optimality criteria used fall into 
two classes: Bayesian and frequent. 


Representing the former approach is the Bayes Risk Criterion. Representing 
the latter is the Neyman-Pearson Criterion. These two approaches are 
developed at length in separate modules. 


Statistics Versus Engineering Lingo 


The following table, adapted from Kay, _p.65, summarizes the different 
terminology for hypothesis testing from statistics and signal processing: 


Statistics Signal Processing 


Hypothesis Test Detector 

Null Hypothesis Noise Only Hypothesis 
Alternate Hypothesis Signal + Noise Hypothesis 
Critical Region Signal Present Decision Region 
Type I Error False Alarm 

‘Type II Error Miss 

Size of Test (a) Probability of False Alarm (Pr) 


Power of Test () Probability of Detection (Pp) 


Criteria in Hypothesis Testing 


The criterion used in the previous section - minimize the average cost of an 
incorrect decision - may seem to be a contrived way of quantifying 
decisions. Well, often it is. For example, the Bayesian decision rule depends 
explicitly on the a priori probabilities; a rational method of assigning values 
to these - either by experiment or through true knowledge of the relative 
likelihood of each model - may be unreasonable. In this section, we develop 
alternative decision rules that try to answer such objections. One essential 
point will emerge from these considerations: the fundamental nature of 
the decision rule does not change with choice of optimization criterion. 
Even criteria remote from error measures can result in the likelihood ratio 
test (see this problem). Such results do not occur often in signal processing 
and underline the likelihood ratio test's significance. 


Maximum Probability of a Correct Decision 


As only one model can describe any given set of data (the models are 
mutually exclusive), the probability of being correct P. for distinguishing 
two models is given by 


P. = Prisay -@ when .@ true] + Pr{say 4% when -4 true] 


We wish to determine the optimum decision region placement Expressing 
the probability correct in terms of the likelihood functions p,,_y, (r), thea 


priori probabilities, and the decision regions, 


Po= fm Sa (art [im aad 


We want to maximize P, by selecting the decision regions $%9 and SXg. The 
probability correct is maximized by associating each value of r with the 
largest term in the expression for P,. Decision region So, for example, is 
defined by the collection of values of r for which the first term is largest. 
As all of the quantities involved are non-negative, the decision rule 
maximizing the probability of a correct decision is 


Note:Given r, choose .@; for which the product 7; p,\_g, (7) is largest. 


Simple manipulations lead to the likelihood ratio test. 


Pr|.m™ (7) 2 0 


Pr|H (r) y, 1 


Note that if the Bayes' costs were chosen so that Ci; = 0 and C;; = C, ( 
i #~ 7), we would have the same threshold as in the previous section. 


To evaluate the quality of the decision rule, we usually compute the 
probability of error P. rather than the probability of being correct. This 
quantity can be expressed in terms of the observations, the likelihood ratio, 
and the sufficient statistic. 

Equation: 


Po = To f Pr|m (r)drimf Pr|.m™ (r)dr 


To f Palm (A)d A+ f Paya (A)d A 
= mofprm(Ndvt+mfrera(NHNar 


When the likelihood ratio is non-monotonic, the first expression is most 
difficult to evaluate. When monotonic, the middle expression proves the 
most difficult. Furthermore, these expressions point out that the likelihood 
ratio and the sufficient statistic can be considered a function of the 
observations 7; hence, they are random variables and have probability 
densities for each model. Another aspect of the resulting probability of error 
is that no other decision rule can yield a lower probability of error. This 
statement is obvious as we minimized the probability of error in deriving 
the likelihood ratio test. The point is that these expressions represent a 
lower bound on performance (as assessed by the probability of error). This 
probability will be non-zero if the conditional densities overlap over some 
range of values of r, such as occurred in the previous example. In this 
region of overlap, the observed values are ambiguous: either model is 


consistent with the observations. Our "optimum" decision rule operates in 
such regions by selecting that model which is most likely (has the highest 
probability) of generating any particular value. 


Neyman-Pearson Criterion 


Situations occur frequently where assigning or measuring the a priori 
probabilities P; is unreasonable. For example, just what is the a priori 
probability of a supernova occurring in any particular region of the sky? We 
clearly need a model evaluation procedure which can function without a 
priori probabilities. This kind of test results when the so-called Neyman- 
Pearson criterion is used to derive the decision rule. The ideas behind and 
decision rules derived with the Neyman-Pearson criterion (Neyman and 
Pearson) will serve us well in sequel; their result is important! 


Using nomenclature from radar, where model .4 represents the presence 
of a target and .G@ its absence, the various types of correct and incorrect 
decisions have the following names (Woodward, pp. 127-129).[ footnote] 


e Detection we say it's there when it is; Pp = Pr(say “%|.@% true) 
e False-alarm we say it's there when it's not; 

Pr = Pr(say “@|-%@ true) 
e Miss we say it's not there when it is; Py; = Pr (say |. true) 


The remaining probability Pr[say -Z | Zp true] has historically been 
left nameless and equals 1 — Pr. We should also note that the detection and 
miss probabilities are related by Pjy = 1 — Pp. As these are conditional 
probabilities, they do not depend on the a priori probabilities and the two 
probabilities Pp and Pp characterize the errors when any decision rule is 
used. 

In hypothesis testing, a false-alarm is known as a type I error and a miss a 
type II error. 


These two probabilities are related to each other in an interesting way. 
Expressing these quantities in terms of the decision regions and the 
likelihood functions, we have 


Pe=f Pela (r)dr 


Po= | Pra (r)dr 


As the region $4; shrinks, both of these probabilities tend toward zero; as 
SR expands to engulf the entire range of observation values, they both tend 
toward unity. This rather direct relationship between Pp and Pr does not 
mean that they equal each other; in most cases, as $4; expands, Pp 
increases more rapidly than Pp (we had better be right more often than we 
are wrong!). However, the "ultimate" situation where a rule is always right 
and never wrong (Pp = 1, Pr = 0) cannot occur when the conditional 
distributions overlap. Thus, to increase the detection probability we must 
also allow the false-alarm probability to increase. This behavior represents 
the fundamental tradeoff in hypothesis testing and detection theory. 


One can attempt to impose a performance criterion that depends only on 
these probabilities with the consequent decision rule not depending on the a 
priori probabilities. The Neyman-Pearson criterion assumes that the false- 
alarm probability is constrained to be less than or equal to a specified value 
a while we attempt to maximize the detection probability Pp. 


VPr, Pr < a: (maxm, {R1, Pp}) 


A subtlety of the succeeding solution is that the underlying probability 
distribution functions may not be continuous, with the result that Pr can 
never equal the constraining value a. Furthermore, an (unlikely) possibility 
is that the optimum value for the false-alarm probability is somewhat less 
than the criterion value. Assume, therefore, that we rephrase the 
optimization problem by requiring that the false-alarm probability equal a 
value a’ that is less than or equal to a. 


This optimization problem can be solved using Lagrange multipliers (see 
Constrained Optimization); we seek to find the decision rule that maximizes 


F = Pp+X(Pr-a’) 


where A is the Lagrange multiplier. This optimization technique amounts to 
finding the decision rule that maximizes F’, then finding the value of the 
multiplier that allows the criterion to be satisfied. As is usual in the 
derivation of optimum decision rules, we maximize these quantities with 
respect to the decision regions. Expressing Pp and Pr in terms of them, we 
have 

Equation: 


F = f ppm (r)dr+a(f Prim (nr) dr—a’) 
= —(Aa’) + f Prim (") +A Pom (tr) dr 


To maximize this quantity with respect to $41, we need only to integrate 
over those regions of r where the integrand is positive. The region 9%; thus 
corresponds to those values of r where pr|.~, (r) > — (A pr|.m (r)) and 
the resulting decision rule is 


My, 
z (-A) 
My 


The ubiquitous likelihood ratio test again appears; it is indeed the 
fundamental quantity in hypothesis testing. Using the logarithm of the 
likelihood ratio or the sufficient statistic, this result can be expressed as 
either 


M 
Cees In(—A) 


or 


MM 


Mr)2y 
Mo 


We have not as yet found a value for the threshold. The false-alarm 
probability can be expressed in terms of the Neyman-Pearson threshold in 
two (useful) ways. 


Equation: 


Jy pram (Nar 


One of these implicit equations must be solved for the threshold by setting 
Pr equal to a’. The selection of which to use is usually based on pragmatic 
considerations: the easiest to compute. From the previous discussion of the 
relationship between the detection and false-alarm probabilities, we find 
that to maximize Pp we must allow a’ to be as large as possible while 
remaining less than a. Thus, we want to find the smallest value of —A 
(note the minus sign) consistent with the constraint. Computation of the 
threshold is problem-dependent, but a solution always exists. 


Example: 

An important application of the likelihood ratio test occurs when r is a 
Gaussian random vector for each model. Suppose the models correspond to 
Gaussian random vectors having different mean values but sharing the 
same identity covariance. 


© Myr ~ N (0,071) 
° Mr ~ N(m,o71) 


Thus, 7 is of dimension L and has statistically independent, equal variance 


components. The vector of means m = (mp...mz_1)” distinguishes the 
two models. The likelihood functions associated this problem are 


-(1/2(3)’) 


@=l1 : (eGry) 


The likelihood ratio A(r) becomes 


oe -(1/2(2)') 
A(r) = Atico © 8 
E-1 ,-(1/2(2) } 


[=0 


This expression for the likelihood ratio is complicated. In the Gaussian 
case (and many others), we use the logarithm the reduce the complexity of 
the likelihood ratio and form a sufficient statistic. 

Equation: 


In(A(r)) = —1/2——— We EWE ap Ly ue 
= > ye OTs = es my? 
The likelihood ratio test then has the much simpler, but equivalent form 
al 
yt es (o 2In(n )) +1725 met 
1=0 Mo 


To focus on the model evaluation aspects of this problem, let's assume 
means be equal to a positive constant: m; = m (0).[footnote] 


Note that all that need be known about the observations {77} is their sum. 
This quantity is the sufficient statistic 198 the Gaussian problem: 

Y(r) = do rand y = 07 In() + =. 
Why did the authors assume that the mean was positive? What would 
happen if it were negative? 
When trying to compute the probability of error or the threshold in the 
Neyman-Pearson criterion, we must find the conditional probability 
density of one of the decision statistics: the likelihood ratio, the log- 
likelihood, or the sufficient statistic. The log-likelihood and the sufficient 


Statistic are quite similar in this problem, but clearly we should use the 
latter. One practical property of the sufficient statistic is that it usually 
simplifies computations. For this Gaussian example, the sufficient statistic 
is a Gaussian random variable under each model. 


© mM: V(r) ~ N (0, Lo’) 
© MM: V(r) ~ N (Lm, Lo’) 


To find the probability of error from [link], we must evaluate the area 
under a Gaussian probability density function. These integrals are 
succinctly expressed in terms of Q(a), which denotes the probability that a 
unit-variance, zero-mean Gaussian random variable exceeds x (see 
Probability and Stochastic Processes). As 1 — Q(x) = Q(—2), the 
probability of error can be written as 


mse acl 


An Engi special case occurs when 779 = 1/2 = 7. In this case, 
-y = -" and the probability of error becomes 


nao 2) 


As Q(-) is a monotonically decreasing function, the probability of error 
vim 


decreases with increasing values of the ratio Lm However, as shown in 


this figure, Q(-) decreases in a nonlinear ete Thus, increasing m by a 
factor of two may decrease the probability of error by a larger or a smaller 
factor; the amount of change depends on the initial value of the ratio. 

To find the threshold for the Neyman-Pearson test from the expressions 
given on [link], we need the area under a Gaussian density. 

Equation: 


a Oe) 


As Q(-) is a monotonic and continuous function, we can now set a’ equal 
to the criterion value @ with the result 


7 = VLeQ *(a) 


where Q~'(-) denotes the inverse function of Q(-). The solution of this 
equation cannot be performed analytically as no closed form expression 
exists for Q(-) (much less its inverse function); the criterion value must be 
found from tables or numerical routines. Because Gaussian problems arise 
frequently, the [link] accompanying table provides numeric values for this 
quantity at the decade points. 


x Oma) 
107+ 1.281 
10°? 2.396 
10°° 3.090 
10-4 3.719 
10~° 4.265 
10° 4.754 


The table displays interesting values for Q~'(-) that can be used to 
determine thresholds in the Neyman-Pearson variant of the likelihood ratio 
test. Note how little the inverse function changes for decade changes in its 
argument; Q(-) is indeed very nonlinear. 


The detection probability is given by 


Pp =Q (Xe an 


The Bayes Risk Criterion in Hypothesis Testing 


The design of a hypothesis test/detector often involves constructing the solution 
to an optimization problem. The optimality criteria used fall into two classes: 
Bayesian and frequent. 


In the Bayesian setup, it is assumed that the a priori probability of each 
hypothesis occuring (7;) is known. A cost C}j is assigned to each possible 
outcome: 


Ci; = Pr|sayH;whenH ;true| 


The optimal test/detector is the one that minimizes the Bayes risk, which is 
defined to be the expected cost of an experiment: 
C= S- Ci; Prisay H;whenH true] 
1,J 


In the event that we have a binary problem, and both hypotheses are simple, the 
decision rule that minimizes the Bayes risk can be constructed explicitly. Let us 
assume that the data is continuous (i.e., it has a density) under each hypothesis: 


Ho 7:a~ fo(x) 
Ay ars baad fi(x) 


Let Ro and R, denote the decision regions corresponding to the optimal test. 
Clearly, the optimal test is specified once we specify Ro and R, = Ro’. 


The Bayes risk may be written 
Equation: 


C = yaar Cin [ ia) dz 
Jf Cootofo(#) + Cormfi(a) da + f Cromofo(x) + Cumfi(a) dx 


Recall that Ro and R; partition the input space: they are disjoint and their union 
is the full input space. Thus, every possible input z belongs to precisely one of 
these regions. In order to minimize the Bayes risk, a measurement x should 
belong to the decision region R; for which the corresponding integrand in the 


preceding equation is smaller. Therefore, the Bayes risk is minimized by 
assigning x to Ro whenever 


ToCoofo(x) + mCoifi(@) < moCiofo(#) + mCufi(x) 


and assigning x to R, whenever this inequality is reversed. The resulting rule 
may be expressed concisely as 


_ file) Bh 19 (Cio — Coo) _ 
~ f(z) ® m(Cu—Cu) 


A(a) 


Here, A(a) is called the likelihood ratio, 7) is called the threshold, and the 
overall decision rule is called the Likelihood Ratio Test (LRT). The expression 
on the right is called a threshold. 


Example: 

An instructor in a course in detection theory wants to determine if a particular 
student studied for his last test. The observed quantity is the student's grade, 
which we denote by 7. Failure may not indicate studiousness: conscientious 
students may fail the test. Define the models as 


e .M@: did not study 
© ;: did study 


The conditional densities of the grade are shown in [link]. 


Conditional densities for the grade 
distributions assuming that a student did 
not study (4) or did (.4) are shown in 
the top row. The lower portion depicts the 
likelihood ratio formed from these 
densities. 


Based on knowledge of student behavior, the instructor assigns a priori 
probabilities of 79 = + and 77, = + The costs C;; are chosen to reflect the 
instructor's sensitivity to student feelings: Co, = 1 = C9 (an erroneous 
decision either way is given the same cost) and C99 = 0 = Cj,. The likelihood 
ratio is plotted in [link] and the threshold value 7, which is computed from the a 
priori probabilities and the costs to be + is indicated. The calculations of this 


comparison can be simplified in an obvious way. 


M 
bog i 
50 %, 3 


or 


The multiplication by the factor of 50 is a simple illustration of the reduction of 
the likelihood ratio to a sufficient statistic. Based on the assigned costs and a 
priori probabilities, the optimum decision rule says the instructor must assume 
that the student did not study if the student's grade is less than 16.7; if greater, 
the student is assumed to have studied despite receiving an abysmally low grade 
such as 20. Note that as the densities given by each model overlap entirely: the 
possibility of making the wrong interpretation always haunts the instructor. 
However, no other procedure will be better! 


A special case of the minimum Bayes risk rule, the minimum probability of error 
rule, is used extensively in practice, and is discussed at length in another module. 


Problems 


Exercise: 


Problem: 


Denote a = Pr{declareH;whenHotrue] and 


6 = PrideclareH,whenH; true]. Express the Bayes risk C in terms of a 
and 8, C;;, and 7;. Argue that the optimal decision rule is not altered by 
setting Coo = C11 = 0. 


Exercise: 
Problem: 


Suppose we observe a such that A(a) = 7. Argue that it doesn't matter 
whether we assign zw to Ro or Ry. 


Minimum Probability of Error Decision Rule 
Consider the binary hypothesis test 
KH: x2 ~ fo(x) 
JG, : ar fi(x) 


Let 7;, denote the a priori probability of hypothesis #3. Suppose our 
decision rule declares ".% is the true model" when x € Ro, and it selects 
JE6, when x € R1, where R; = Ro’. The probability of making an error, 
denoted P., is 

Equation: 


P. = Pr{declareAand% true] + Prideclare Mand “true| 
= Pr[ #4] Pr| Ho | 44] + Pr[ Ho] Prl%4 | H] 
f mfi(x)da+ f mofo(x) dx 


| 


In this module, we study the minimum probability of error decision rule, 
which selects Rg and R; so as to minimize the above expression. 


Since an observation @ falls into one and only one of the decision regions 
R;, in order to minimize P., we assign z to the region for which the 
corresponding integrand in [link] is smaller. Thus, we select x € Ro if 
™fi(x) < mofo(a), and x € R, if the inequality is reversed. This 
decision rule may be summarized concisely as 


— Ale) J mm _ 


~ fo(e) 3% ™ 


A(z) 


Here, A(a) is called the likelihood ratio, 7 is called a threshold, and the 
overall decision rule is called the Likelihood Ratio Test. 


Example: 


Normal with Common Variance, Uncommon 
Means 


Consider the binary hypothesis test of a scalar x 
Pe ee (0, a”) 
TON (u, a”) 
where p and o” are known, positive quantities. Suppose we 


observe a single measurement zx. The likelihood ratio is 
Equation: 


Ae) = = 


The expression for A(x) is somewhat complicated. By 
applying a sequence of monotonically increasing functions 
to both sides, we can obtain a simplified expression for the 
optimal decision rule without changing the rule. In this 
example, we apply the natural logarithm and rearrange 
terms to atrive at 


= 


Here we have used the assumption pz > 0. If  < 0, then 
dividing by yz would reverse the inequalities. 


This form of the decision rule is much simpler: we just 
compare the observed value z to a threshold 7. [link] 
depicts the two candidate densities and a possible value of y 
. If each hypothesis is a priori equally likely (7) = 71 = 
), then y = 5. [link] illustrates the case where 779 > 771 ( 
y> $). 


The two candidate densities, and a threshold 
corresponding to 7 > 7 


If we plot the two densities so that each is weighted by its a 
priori probability of occuring, the two curves will intersect 
at the threshold +¥ (see [link]). (Can you explain why this is? 
Think back to our derivation of the LRT). This plot also 
offers a way to visualize the probability of error. Recall 
Equation: 


P. = fmfi(z)dz2+ f mofo(z) da 
( eagenee) clas 5 | aacin(e) clas 
= Wily + torr 


where Py and Pr denote the miss and false alarm 
probabilities, respectively. These quantities are depicted in 
[link]. 


The candidate densities weighted by their a priori 
probabilities. The shaded region is the probability of 
error for the optimal decision rule. 


We can express Py and Pp in terms of the Q-function as 


r=no(#22) +me(2) 


When 779 = 71 = ot we have 7 = . , and the error 


probability is 


r= a(t) 


Since Q(a) is monotonically decreasing, this says that the 
"difficulty" of the detection problem decreases with 
decreasing o and increasing ju. 


In the preceding example, computation of the probability of error involved 
a one-dimensional integral. If we had multiple observations, or vector- 
valued data, generalizing this procedure would involve multi-dimensional 
integrals over potentially complicated decision regions. Fortunately, in 
many cases, we can avoid this problem through the use of sufficient 
Statistics. 


Example: 
Suppose we have the same test as in the previous example, but now we 
have NV independent observations: 


Foy oy (Oo \h eee 
Ae Oey (tao |e lar 


where ps > 0 and o? > 0 and both are known. The likelihood ratio is 
Equation: 


= (an—p)? 


N 1 
= e 20 
A(a) —, Ilha V2n02 
— ao 
a ee 
= V 2a 


ss 

e 202 ee (en—p)? 

“1 5N .,,2 
edt ~n=1" 

1 N 2 

= 2 Lin=1 2%np— 


1 N Np2 
— .fle tees 


As in the previous example, we may apply the natural logarithm and 
rearrange terms to obtain an equivalent form of the LRT: 


The scalar quantity ¢ is a sufficient statistic for the mean. In order to 
evaluate the probability of error without resorting to a multi-dimensional 
integral, we can express P, in terms of ¢ as 


P. = 7 Prit <y| Mtrue] + 79 Prit > y | Mtruel 


Now t is a linear combination of normal variates, so it is itself normal. In 
particular, we have t = Aw, where (1 ... 1) is an N-dimensional row 


vector of 1's, and & is multivariate normal with mean O or pp = (p...)’, 
and covariance o7J. Thus we have 


t| 4 ~ WN (AO, Ao? IA‘) = N (0, No?) 
t| 4, ~ N (Ap, Ao*IA*) = W(Nu, No’) 


Therefore, we may write P. in terms of the Q-function as 


P= m(AB—7) + m9( ul ) 
No No 


In the special case 7 = 7 = 5, 


Since Q is monotonically decreasing, this result provides mathematical 
support for something that is intuitively obvious: The performance of our 
decision rule improves with increasing NV and p, and decreasing o. 


Note:In the context of signal processing, the foregoing problem may be 
viewed as the problem of detecting a constant (DC) signal in additive 


white Gaussian noise: 
Ty Lie Watt — la 
TOL, = A+ Wy = Lang VN 


where A is a known, fixed amplitude, and w, ~ VW (0, Gai Here A 
corresponds to the mean yz in the example. 


The next example explores the minimum probability of error decision rule 
in a discrete setting. 


Example: 
Repetition Code 


Suppose we have a friend who is trying to transmit a bit (0 
or 1) to us over a noisy channel. The channel causes an 
error in the transmission (that is, the bit is flipped) with 
probability p, where 0 < p < + and p is known. In order 
to increase the chance of a successful transmission, our 
friend sends the same bit NV times. Assume the VV 
transmissions are statistically independent. Under these 
assumptions, the bits you receive are Bernoulli random 
variables: 2, ~ Bernoulli (0). We are faced with the 
following hypothesis test: 


. A—m 1’ cant 


uly vu —_ vp VU OCLit 


KH d=1-p 1 sent 


We decide to decode the received sequence 
Deon nv)" by minimizing the probability of error. 
The likelihood ratio is 

Equation: 


tn, 1-xr2n 


a ESS ee 
Ae) [[n sp?" (1—p)* 
(1—p)*p-* 


where k = San Ly, is the number of 1s received. 
Note:k is a sufficient statistic for 0. 


The LRT is 


low 2k-N a 
Cs) E 
Pp Ha M1 
Taking the natural logarithm of both sides and rearranging, 
we have 


N 1 In 

N 1 _In(n)_ 

2 2 In( +2) 
Pp 


In the case that both hypotheses are equally likely, the 


H, 
k2 
Hy 


minimum probability of error decision is the "majority- 
vote" rule: Declare .% if there are more 1s than Os, declare 
Ay otherwise. In the event k = -y, we may decide 
arbitrarily; the probability of error is the same either way. 
Let's adopt the convention that #% is declared in this case. 


To compute the probability of error of the optimal rule, 
write 
Equation: 


P, = moPrideclareH | Motrue] + 7 Prideclare% | Mtr 
= moPrik >>| #otrue] + 7 Pr[k < y | Mtrue] 


Now k is a binomial random variable, 

k, ~ Binomial (NV, 8), where @ depends on which 
hypothesis is true. We have 

Equation: 


Pri eral Se fo(k) 


N = 
= ees e a —p)"* 
and 


ly] 
Prik <7 #4] = 90 (), )(a— yp 


k=0 


Using these formulae, we may compute P, explicitly for 
given values of N, p, 7 and 7. 


MAP Interpretation 


The likelihood ratio test is one way of expressing the minimum probability 
of error decision rule. Another way is 
Rule 


Declare hypothesis 7 such that 7; f;(a) is maximal. 

This rule is referred to as the maximum a posteriori, or MAP rule, 
because the quantity 7; f;(a) is proportional to the posterior probability of 
hypothesis 7. This becomes clear when we write 7; = Pr|.#%4] and 

fi(x) = f(a|#). Then, by Bayes rule, the posterior probability of 
given the data is 


Pr| 41] f(2| Hi) 
f(z) 


Here f(a) is the unconditional density or mass function for a, which is 
effectively a constant when trying to maximiaze with respect to 2. 


Pr[3G | @] = 


According to the MAP interpretation, the optimal decision boundary is the 
locus of points where the weighted densities (in the continuous case) 
1; f;() intersect one another. This idea is illustrated in [link]. 


Multiple Hypotheses 


One advantage the MAP formulation of the minimum probability of error 
decision rule has over the LRT is that it generalizes easily to M/-ary 
hypothesis testing. If we are to choose between hypotheses 7%, 

i = {1,..., M}, the optimal rule is still the MAP rule 


Special Case of Bayes Risk 


The Bayes risk criterion for constructing decision rules assigns a cost C’,; to 
the outcome of declaring .# when .#; is in effect. The probability of error 
is simply a special case of the Bayes risk corresponding to Cog = Cy, = 0 
and Co, = Cio = 1. Therefore, the form of the minimum probability of 


error decision rule is a specialization of the minimum Bayes risk decision 
rule: both are likelihood ratio tests. The different costs in the Bayes risk 
formulation simply shift the threshold to favor one hypothesis over the 
other. 


Problems 


Exercise: 
Problem: 
Generally speaking, when is the probability of error zero for the 
optimal rule? Phrase your answer in terms of the distributions 


underlying each hypothesis. Does the LRT agree with your answer in 
this case? 


Exercise: 
Problem: 
Suppose we measure NV independent values 71,..., “1. We know the 


variance of our measurements (a? = 1), but are unsure whether the 
data obeys a Laplacian or Gaussian probability law: 


Hs: fola) = a 


1 2 
Hr tie) = aa 
TT 


Show that the two densities have the same mean and variance, and plot 
the densities on the same graph. 


Find the likelihood ratio. 


Determine the decision regions for different values of the threshold 7. 
Consider all possible values of 7 > 0 


Note:There are three distinct cases. 


Draw the decision regions and decision boundaries for 7 = oe 1G 2\. 


Assuming the two hypotheses are equally likely, compute the 
probability of error. Your answer should be a number. 


Exercise: 


Problem: 


Arbitrary Means and Covariances 


Consider the hypothesis testing problem 
KH a 67 ~ NV (pu, Xo) 
ia 7 & ~ MN (i, 1) 


where jo € IR? and Lo € R¢, and 2/9, 51 are positive definite, 
symmetric dxd matrices. Write down the likelihood ratio test, and 
simplify, for the following cases. In each case, provide a geometric 
description of the decision boundary. 


dig = 41, but lo F H1. 


Lo = pi, but Xo A DY. 


Mo A fy and Lig F LX. 
Exercise: 


Problem: 


Suppose we observe WN independent realizations of a Poisson random 
variable k with intensity parameter A: 


e>yk 


f(k) = 


We must decide which of two intensities is in effect: 
FO XS 
3: rX=X1 

where Ap < Aj. 

Give the minimum probability of error decision rule. 


Simplify the LRT to a test statistic involving only a sufficient statistic. 
Apply a monotonically increasing transformation to simplify further. 


Determine the distribution of the sufficient statistic under both 
hypotheses. 


Note: Use the characteristic function to show that a sum of IID 
Poisson variates is again Poisson distributed. 


Derive an expression for the probability of error. 


Assuming the two hypotheses are equally likely, and Ag = 5 and 
A 1 = 6, what is the minimum number NV of observations needed to 
attain a probability of error no greater than 0.01? 


Note:If you have numerical trouble, try rewriting the log-factorial so 
as to avoid evaluating the factorial of large integers. 


Exercise: 


Problem: 


In [link], suppose 7) = 7, = + and p = 0.1. What is the smallest 
value of NV needed to ensure P. < 0.01? 


The Neyman-Pearson Criterion 


In hypothesis testing, as in all other areas of statistical inference, there are two 
major schools of thought on designing good tests: Bayesian and frequentist (or 
classical). Consider the simple binary hypothesis testing problem 


KH: 2 ~ fo(x) 
FG: a~ fi(x) 


In the Bayesian setup, the prior probability 7; = Pr|.#] of each hypothesis 
occurring is assumed known. This approach to hypothesis testing is 
represented by the minimum Bayes risk criterion and the minimum probability 
of error criterion. 


In some applications, however, it may not be reasonable to assign an a priori 
probability to a hypothesis. For example, what is the a priori probability of a 
supernova occurring in any particular region of the sky? What is the prior 
probability of being attacked by a ballistic missile? In such cases we need a 
decision rule that does not depend on making assumptions about the a priori 
probability of each hypothesis. Here the Neyman-Pearson criterion offers an 
alternative to the Bayesian framework. 


The Neyman-Pearson criterion is stated in terms of certain probabilities 
associated with a particular hypothesis test. The relevant quantities are 
summarized in [link]. Depending on the setting, different terminology is used. 


Statistics Signal Processing 

Probability Name Notation Name Notation 
false- 

Po(declareH ) size a alarm Pr 


probability 


Statistics Signal Processing 


Probability Name Notation Name Notation 
detection 
P, (declare ) power B arobability Pp 


Here P;(declare.#%;) dentoes the probability that we declare hypothesis #4 to 
be in effect when #% is actually in effect. The probabilities Po (declare) 
and P; (declare) (sometimes called the miss probability), are equal to 

1 — Pr and 1 — Pp, respectively. Thus, P- and Pp represent the two degrees 
of freedom in a binary hypothesis test. Note that Pp and Pp do not involve a 
priori probabilities of the hypotheses. 


These two probabilities are related to each other through the decision regions. 
If R, is the decision region for #4, we have 


Pe= f folw)de 


Pp =f Aw dz 


The densities f;(a) are nonnegative, so as R, shrinks, both probabilities tend 
to zero. As R; expands, both tend to one. The ideal case, where Pp = 1 and 
Py = 0, cannot occur unless the distributions do not overlap (i.e., 

f fo(x) fi(a) d x = 0). Therefore, in order to increase Pp, we must also 
increase Pp. This represents the fundamental tradeoff in hypothesis testing 
and detection theory. 


Example: 
Consider the simple binary hypothesis test of a scalar measurement x: 


Foe ey (0,1) 
oi ee (lee) 


Suppose we use a threshold test 


where 7 € R is a free parameter. Then the false alarm and detection 


probabilities are 
CO 
1 
Pe= [ 
y V2T 


sae | (20E 
Pp= | e °F? di=Q(y-1) 
y V/27 


where @ denotes the Q-function. These quantities are depicted in [link]. 


eo? dt =Q(1) 


False alarm and detection values for a certain threshold. 


Since the Q-function is monotonicaly decreasing, it is evident that both Pp 
and Pr decay to zero as ¥y increases. There is also an explicit relationship 


Pp = Q(Q" (Pr) — 1) 


A common means of displaying this relationship is with a receiver operating 
characteristic (ROC) curve, which is nothing more than a plot of Pp versus 
Pr ((ink)). 


ROC curve for this example. 


The Neyman-Pearson Lemma: A First Look 


The Neyman-Pearson criterion says that we should construct our decision rule 
to have maximum probability of detection while not allowing the probability 
of false alarm to exceed a certain value a. In other words, the optimal detector 


according to the Neyman-Pearson criterion is the solution to the following 


constrainted optimization problem: 


Neyman-Pearson Criterion 


Equation: 


max {Pp},such thatPr < a 


The maximization is over all decision rules (equivalently, over all decision 
regions Ro, R,). Using different terminology, the Neyman-Pearson criterion 
selects the most powerful test of size (not exceeding) a. 


Fortunately, the above optimization problem has an explicit solution. This is 
given by the celebrated Neyman-Pearson lemma, which we now state. To 
ease the exposition, our initial statement of this result only applies to 
continuous random variables, and places a technical condition on the densities. 
A more general statement is given later in the module. 

Neyman-Pearson Lemma: initial statement 


Consider the test 


fi(z) 
fo(x) 
satisfies the condition that for each y € R, A(a) takes on the value y with 
probability zero under hypothesis .. The solution to the optimization 
problem in [link] is given by 


where f;(a) is a density. Define A(a) = , and assume that A(a:) 


where 77 is such that 
Pr= [ foe) de=a 


If ~ = 0, then 7 = oo. The optimal test is unique up to a set of probability 
zero under #% and “%. 

The optimal decision rule is called the likelihood ratio test. A(x) is the 
likelihood ratio, and 7 is a threshold. Observe that neither the likelihood 


ratio nor the threshold depends on the a priori probabilities Pr[.#]. they 
depend only on the conditional densities f; and the size constraint a. The 
threshold can often be solved for as a function of a, as the next example 
shows. 


Example: 
Continuing with [link], suppose we wish to design a Neyman-Pearson 
decision rule with size constraint a. We have 


Equation: 
eee eat 
Vir 
A(z) = = 
a 
— pio 


By taking the natural logarithm of both sides of the LRT and rarranging 
terms, the decision rule is not changed, and we obtain 


FE, 1 
z2In(n)+ 5 = 
KH 


Thus, the optimal rule is in fact a thresholding rule like we considered in 
[link]. The false-alarm probability was seen to be 


Pr = Q(y) 


Thus, we may express the value of -y required by the Neyman-Pearson lemma 
in terms of a: 


y= Q" (a) 


Sufficient Statistics and Monotonic Transformations 


For hypothesis testing involving multiple or vector-valued data, direct 
evaluation of the size (Pr) and power (Pp) of a Neyman-Pearson decision 
rule would require integration over multi-dimensional, and potentially 
complicated decision regions. In many cases, however, this can be avoided by 
simplifying the LRT to a test of the form 


where the test statistic £ = T(a:) is a sufficient statistic for the data. Such a 
simplified form is arrived at by modifying both sides of the LRT with 
montonically increasing transformations, and by algebraic simplifications. 
Since the modifications do not change the decision rule, we may calculate Pp 
and Pp in terms of the sufficient statistic. For example, the false-alarm 
probability may be written 

Equation: 


Pr = Prideclare#4| 
= f fo(t)dt 


where f(t) denotes the density of ¢ under @. Since t is typically of lower 
dimension than 2, evaluation of Pr and Pp can be greatly simplified. The key 
is being able to reduce the LRT to a threshold test involving a sufficient 
statistic for which we know the distribution. 


Example: 
Common Variances, Uncommon Means 


Let's design a Neyman-Pearson decision rule of size a for the 
problem 


KH, : a~ NV (0,071) 
UG eae ~ NV (pl, o71) 


where pz > 0, 0? > O are known, 0 = (0.. ‘Oe 


1 = (1...1)* are N-dimensional vectors, and I is the Nx N 
identity matrix. The likelihood ratio is 


Equation: 
fe , _ (en-1)? 
Qo2 
A a Tn Ino : 
(x) = 2 
N 1 Ss 
Tn V2n02 a 
(tn—")? 
= ene Ee 
ro 2 
e one oe 


N 
— ae = ast an p— pe 
N, N 
+5 (- £ +p eee zy) 


To simplify the test further we may apply the natural logarithm 
and rearrange terms to obtain 


Note: We have used the assumption p > 0. If ~ < 0, then 
division by p is not a monotonically increasing operation, and 
the inequalities would be reversed. 


The test statistic t is sufficient for the unknown mean. To set 
the threshold y, we write the false-alarm probability (size) as 


a= eae ol = / fo(t) dt 


To evaluate Pr, we need to know the density of ¢ under %. 
Fortunately, t is the sum of normal variates, so it is again 


normally distributed. In particular, we have t = Aa, where 
A=1' so 


t ~ W(AO, A (o*I)A*) = V(0, No”) 


under .#%. Therefore, we may write Pr in terms of the Q- 
function as 


ro He) 


The threshold is thus determined by 
y= VNoQ" (a) 
Under #4, we have 
t~ W(Al1, A (o7I)A*) = WY (Ny, No’) 


and so the detection probability (power) is 


Writing Pp as a function of Pr, the ROC curve is given by 


2p = o(orrn = 2M) 


The quantity iE is called the signal-to-noise ratio. As its 
name suggests, a larger SNR corresponds to improved 


performance of the Neyman-Pearson decision rule. 


Note:In the context of signal processing, the foregoing 
problem may be viewed as the problem of detecting a 
constant (DC) signal in additive white Gaussian noise: 


Go) a SN, 
(ses peated as 11g etal Deen A 


where A is a known, fixed amplitude, and w, ~ (0,07). 
Here A corresponds to the mean yp in the example. 


The Neyman-Pearson Lemma: General Case 


In our initial statement of the Neyman-Pearson Lemma, we assumed that for 
all n, the set {a, w| A(a) = 7} had probability zero under #. This 
eliminated many important problems from consideration, including tests of 
discrete data. In this section we remove this restriction. 


It is helpful to introduce a more general way of writing decision rules. Let y 
be a function of the data # with y(a) € [0, 1]. y defines the decision rule 
"declare #4, with probability y(a)." In other words, upon observing az, we 
flip a "p(a) coin." If it turns up heads, we declare #4; otherwise we declare 
4. Thus far, we have only considered rules with y(a) € {0, 1} 
Neyman-Pearson Lemma 


Consider the hypothesis testing problem 
KH: 2 ~ fo(a) 
A: a ~ fi(x) 


where fg and f; are both pdfs or both pmfs. Let a € [0, 1) be the size (false- 
alarm probability) constraint. The decision rule 


is the most powerful test of size a, where 7 and p are uniquely determined by 
requiring Pp = a. If a = 0, we take 7 = oo, p = 0. This test is unique up to 
sets of probability zero under #% and #%. 

When Pr|A(a) = 7] > 0 for certain 7, we choose 7 and p as follows: Write 


Pp = Pr[A(x) > yn] + p Pr[ A(x) = 7] 


Choose 77 such that 


Then choose p such that 


p Pr|A(x) = 7] = a — Pr[A(@) < 7] 


Example: 
Repetition Code 


Suppose we have a friend who is trying to transmit a bit (0 or 
1) to us over a noisy channel. The channel causes an error in 
the transmission (that is, the bit is flipped) with probability p, 
where 0 < p < 4, and pis known. In order to increase the 
chance of a successful transmission, our friend sends the same 
bit NV times. Assume the NV transmissions are statistically 
independent. Under these assumptions, the bits you receive are 
Bernoulli random variables: x, ~ Bernoulli (0). We are faced 
with the following hypothesis test: 


KH : 0 = p(0 sent) 
JE, : 0 = 1 -— p(1 sent) 


We decide to decode the received sequence x = (#1...2 ny) 
by designing a Neyman-Pearson rule. The likelihood ratio is 
Equation: 


Tin. =p)*"pl-*" 
Ae2) TIN, p?»(1—p)*" 
(1—p)*p¥-* 


where k = )~"_, x, is the number of 1s received. 


Note:k is a sufficient statistic for 0. 


The LRT is 


Taking the natural logarithm of both sides and rearranging, we 
have 

SN 1 Inf) 
Se aye 


k 


The false alarm probability is 
Equation: 


Pr = Prlk> | + pPrlk= | 
a an (7 Joka =p o( \ora =H! 


y and p are chosen so that Pr = a, as described above. 


The corresponding detection probability is 


Equation: 


Ey = lee a Salers =) 


Ebest (:) Oe C (1—p)"p" 


Problems 


Exercise: 


Problem: 


Design a hypothesis testing problem involving continous random 
variables such that Pr[A(x) = 7] > 0 for certain values of 7. Write down 
the false-alarm probability as a function of the threshold. Make as general 
a statement as possible about when the technical condition is satisfied. 


Exercise: 
Consider the scalar hypothesis testing problem 
KH: 2 ~ fo(x) 
KA 2 ~ fi(z) 
where 
Problem: f(x) = a a =40, 1} 


n (1+ (2-7) 


Write down the likelihood ratio test. 


Determine the decision regions as a function of 7, for all 7 > 0. Draw a 
representative of each. What are the "critical" values of 7? 


Note:There are five distinct cases. 


Compute the size and power (Pr and Pp) in terms of the threshold 7); 
and plot the ROC. 


Note: 


1 
———_—_ (7 — arctan 2 
eer @) 


Ral 
Suppose we decide to use a simple threshold test x 2 7 instead of the 
KH 


Neyman-Pearson rule. Does our performance # suffer much? Plot the 
ROC for this decision rule on the same graph as for the previous ROC. 


Exercise: 


Problem: 


Suppose we observe NV independent realizations of a Poisson random 
variable k with intensity parameter A: 
ee 

k! 


We must decide which of two intensities is in effect: 
KH A= Ao 
al A= AI 


where Ag < Aj. 


Write down the likelihood ratio test. 


Simplify the LRT to a test statistic involving only a sufficient statistic. 
Apply a monotonically increasing transformation to simplify further. 


Determine the distribution of the sufficient statistic under both 
hypotheses. 


Note: Use the characteristic function to show that a sum of IID Poisson 
variates is again Poisson distributed. 


Derive an expression for the probability of error. 


Assuming the two hypotheses are equally likely, and Ag = 5 and A; = 6, 
what is the minimum number WV of observations needed to attain a false- 
alarm probability no greater than 0.01? 


Note:If you have numerical trouble, try rewriting the log-factorial so as 
to avoid evaluating the factorial of large integers. 


Exercise: 


Problem: 


In [link], suppose p = 0.1. What is the smallest value of NV needed to 
ensure Pr < 0.01? What is Pp in this case? 


The Minimum Variance Unbiased Estimator 


In Search of a Useful Criterion 


In parameter estimation, we observe an N-dimensional vector X of 
measurements. The distribution of X is governed by a density or 
probability mass function fg (a), which is parameterized by an unknown 
parameter @. We would like to establish a useful criterion for guiding the 


_—e 
design and assessing the quality of an estimator 0(a). We will adopt a 
classical (frequentist) view of the unknown parameter: it is not itself 
random, it is simply unknown. 


One possibility is to try to design an estimator that minimizes the mean- 
squared error, that is, the expected squared deviation of the estimated 
parameter value from the true parameter value. For a scalar parameter, the 
MSE is defined by 

Equation: 


For a vector parameter @, this definition is generalized by 
Equation: 


MSE (6,0) = B| (5e) - 6) (0(2) - 0)| 


The expectation is with respect to the distribution of X. Note that for a 
given estimator, the MSE is a function of 0. 


While the MSE is a perfectly reasonable way to assess the quality of an 
estimator, it does not lead to a useful design criterion. Indeed, the estimator 
that minimizes the MSE is simply the estimator 

Equation: 


Unfortunately, this depends on the value of the unknown parameter, and is 
therefore not realizeable! We need a criterion that leads to a realizeable 
estimator. 


Note: In the Bayesian Approach to Parameter Estimation, the MSE is a 
useful design rule. 


The Bias-Variance Decomposition of the MSE 


It is possible to rewrite the MSE in such a way that a useful optimality 
criterion for estimation emerges. For a scalar parameter 9, [Insert 1] This 
expression is called the bias-variance decomposition of the mean-squared 
error. The first term on the right-hand side is called the variance of the 
estimator, and the second term on the right-hand side is the square of the 
bias of the estimator. The formal definition of these concepts for vector 
parameters is now given: 


Let @ be an estimator of the parameter 0. 


variance 
The variance of @ is [Insert 2] 


bias 
The bias of @ is [Insert 3] 


The bias-variance decomposition also holds for vector parameters: [Insert 
4] The proof is a straighforward generalization of the argument for the 
scalar parameter case. 

Exercise: 


Problem: 


Prove the bias-variance decomposition of the MSE for the vector 
parameter case. 


The Bias-Variance Tradeoff 


The MSE decomposes into the sum of two non-negative terms, the squared 
bias and the variance. In general, for an arbitrary estimator, both of these 
terms will be nonzero. Furthermore, as an estimator is modified so that one 
term increases, typically the other term will decrease. This is the so-called 
bias-variance tradeoff. The following example illustrates this effect. 


Example: 


eA an er Ln, Where 2, = A+ Wn, Wn ~ AAO, bi and a is 
an arbitrary constant. 

Let's find the value of a that minimizes the MSE. 

Equation: 


Note: A = aS, Sy ~ N (A, $ | 


Equation: 


- 9 , 

MSE (4) = B| A | —2B|A]A+ A? 
a a7 B| te eee | — 2aB| + ora in| A+ A? 
_ vt a E|x;x;| — 2a x + Se ot Elz,] + A? 

A? +o? if i=j 
A= god 
ME at 
Equation: 
MSE(A) = a?(A?+%) —204?+ 4? 
= ae 4s (a = 1)? A? 
3D aza2 
o(A) =e 
Bias” (4) =(o= ne 
&@ MSE(A ; 
2 
( Se Ig A 0 
Oa 
Equation: 
+ A’ 
) Snes 
eons 


The optimal value a” dpends on the unknown parameter A! Therefore the 
estimator is not realizable. 

Note that the problematic dependence on the parameter enters through the 
Bias component of the MSE. Therefore, a reasonable alternative is to 
constrain the estimator to be unbiased, and then find the estimator that 


produces the minimum variance (and hence provides the minimum MSE 
among all unbiased estimators). 


Note:Sometimes no unbiased estimator exists, and we cannot proceed at 
all in this direction. 


In this example, note that as the value of a@ varies, one of the squared bias 
or variance terms increases, while the other one decreases. Futhermore, note 
that the dependence of the MSE on the unknown parameter is 
manifested in the bias. 


Unbiased Estimators 


Since the bias depends on the value of the unknown parameter, it seems that 
any estimation criterion that depends on the bias would lead to an 
unrealizable estimator, as the previous example suggests (although in 
certain cases realizable minimum MSE estimators can be found). As an 
alternative to minimizing the MSE, we could focus on estimators that have 
a bias of zero. In this case, the bias contributes zero to the MSE, and in 
particular, it does not involve the unknown parameter. By focusing on 
estimators with zero bias, we may hope to arrive at a design criterion that 
yields realizable estimators. 


unbiased 


An estimator 0 is called unbiased if its bias is zero for all values of the 
unknown parameter. Equivalently, [Insert 5] 


For an estimator to be unbiased we require that on average the estimator 
will yield the true value of the unknown parameter. We now give some 
examples. 


The sample mean of a random sample is always an unbiased estimator for 
the mean. 


Example: 
Estimate the DC level in the Guassian white noise. 
Suppose we have data %1,..., Z~ and model the data by 


Vai ete (fq — Ao 
where A is the unknown DC level, and w,, ~ VY (c, a) 


The parameter is —co < A < ow. 
Consider the sample-mean estimator: 


ee ‘le 
Boe 


Is A unbiased? Yes. 
Since E|-] is a linear operator, 


Therefore, A is unbiased! 


What does the unbiased restriction really imply? Recall that g= g(x), a 
function of the data. Therefore, 


and 


Hence, to be unbiased, the estimator (g(-)) must satisfy an integral 
equation involving the densities p (x|8). 
It is possible that an estimator can be unbiased for some parameter values, 


but be biased for others. 


The bias of an estimator may be zero for some values of the unknown 
parameter, but not others. In this case, the estimator is not an unbiased 


estimator. 


Example: 
, lime 
A= — 3 
2N ae 


i _1,_ [0 if (A=0) = unbiased 
a > 5A if (A#0) = biased 


An unbiased estimator is not necessarily a good estimator. 
Some unbiased estimators are more useful than others. 


Example: 


VYWn, Wn ~ No, o*) a(t =A aw) 


= 

B|A,| =A 
o(4,) =o? 
ay'-2 


Both estimators are unbiased, but A» has a much lower variance and 
therefore is a better estimator. 


—_— aS, 


Note: A; (JV) is inconsistent. A2(JV) is consistent. 


Minimum Variance Unbiased Estimators 


Direct minimization of the MSE generally leads to non-realizable 
estimators. Since the dependence of an estimator on the unknown parameter 
appears to come from the bias term, we hope that constraining the bias to be 
zero Will lead to a useful design criterion. But if the bias is zero, then the 
mean-squared error is just the variance. This gives rise to the minimum 
variance unbiased estimator (MVUE) for 8. 


MVUE 


An estimator @ is the minimum variance unbiased estimator if it is 
unbiased and has the smallest variance of any unbiased estimator for 


all values of the unknown parameter. In other words, the MVUE 
satisfies the following two properties: [Insert 6] 


The minimum variance unbiased criterion is the primary estimation 
criterion in the classical (non-Bayesian) approach to parameter estimation. 
Before delving into ways of finding the MVUE, let's first consider whether 
the MVUE always exists. 


Existence of the MVUE 


The MVUE does not always exist. In fact, it may be that no unbiased 
estimators exist, as the following example demonstrates. 


Place [Insert 7] here and make it an example (5). 


Even if unbiased estimators exist, it may be that no single unbiased 
estimator has the minimum variance for all values of the unknown 
parameter. 


Place [Insert 8] here and make it an example (6). 
Exercise: 


Problem: 


Compute the variances of the estimators in the previous examples. 
Using the Cramer-Rao Lower bound, show that one of these two 
estimators has minimum variance among all unbiased estimators. 
Deduce that no single realizable estimator can have minimum variance 
among all unbiased estimators for all parameter values (i.e., the 
MVUE does not exist). When using the Cramer-Rao bound, note that 
the likelihood is not differentable at 9 = 0. 


Methods for Finding the MVUE 


Despite the fact that the MVUE doesn't always exist, in many cases of 
interest it does exist, and we need methods for finding it. Unfortunately, 
there is no 'turn the crank' algorithm for finding MVUE's. There are, 


instead, a variety of techniques that can sometimes be applied to find the 
MVUE. These methods include: 


1. Compute the Cramer-Rao Lower Bound, and check the condition for 
equality. 

2. Find a complete sufficient statistic and apply the Rao-Blackwell 
Theorem. 

3. If the data obeys a general linear model, restrict to the class of linear 
unbiased estimators, and find the minimum variance estimator within 
that class. This method is in general suboptimal, although when the 
noise is Gaussian, it produces the MVUE. 


The Cramer-Rao Lower Bound 


The Cramer-Rao Lower Bound (CRLB) sets a lower bound on the variance of any unbiased 
estimator. This can be extremely useful in several ways: 


1. If we find an estimator that achieves the CRLB, then we know that we have found an MVUB 
estimator! 

2. The CRLB can provide a benchmark against which we can compare the performance of any 
unbiased estimator. (We know we're doing very well if our estimator is "close" to the CRLB.) 

3. The CRLB enables us to rule-out impossible estimators. That is, we know that it is physically 
impossible to find an unbiased estimator that beats the CRLB. This is useful in feasibility 
studies. 

4. The theory behind the CRLB can tell us if an estimator exists that achieves the bound. 


Estimator Accuracy 


Consider the likelihood function p (x|0), where 0 is a scalar unknown (parameter). We can plot the 
likelihood as a function of the unknown, as shown in [link]. 


[missing resource: | 


The more "peaky" or "spiky" the likelihood function, the easier it is to determind the unknown 
parameter. 


Example: 
Suppose we observe 


z£=A+w 


where w ~ (0,07) and A is an unknown parameter. The "smaller" the noise w is, the easier it 
will be to estimate A from the observation x. 

Suppose A = 3ando = 1/3. 

[missing_resource: ] 

Given this density function, we can easily rule-out estimates of A greater than 4 or less than 2, 
since it is very unlikely that such A could give rise to out observation. 

On the other hand, suppose o = 1. 

[missing_resource: ] 

In this case, it is very difficult to estimate A. Since the noise power is larger, it is very difficult to 
distinguish A from the noise. 

The key thing to notice is that the estimation accuracy of A depends on a, which in effect 
determines the peakiness of the likelihood. The more peaky, the better localized the data is about 
the true parameter. 

To quantify the notion, note that the peakiness is effectively measured by the negative of the second 
derivative of the log-likelihood at its peak, as seen in [link]. 

[missing_resource: ] 


Example: 


z=A+w 
Equation: 
1 
log p (x| A) = (- log V 2no?) Sat = )? 
20 
Ol A 
og p («| A) ee ees 
OA oe 
Equation: 


O7log p (z| A) 1 


Omsup o? 


The curvature increases as a” decreases (curvature=peakiness). 


O*logp(z | A) 
Omsup 


In general, the curavture will depend on the observation data; — is a function of x. 


Therefore, an average measure of curvature is more appropriate. 
Equation: 


This average-out randomness due to the data and is a function of @ alone. 


We are now ready to state the CRLB theorem. 
Cramer-Rao Lower Bound Theorem 


Assume that the pdf p (| 6) satisfies the "regularity" condition 


(of 


where the expectation is take with respect to p (a |). Then, the variance of any unbiased estimator 


6 must satisfy 
Equation: 


a\ 2 1 
o(4) 2 —E| yee!) 


where the derivative is evaluated at the true value of 6 and the expectation is with respect to p (x|6) 


. Moreover, an unbiased estimator may be found that attains the bound for all 0 if and only if 
Equation: 


Dlog p (z|) _ 
——ag = 29) (918) — 8) 


for some functions g and I. 


The corresponding estimator is MVUB and is given by 6= g(x), and the minimum variance is 


7(8) 


Example: 


x=A-+w 


where 


oe e\ B a AAD 
Therefore, any unbiased estimator A has o (4) > o”. But we know that A = z has o (4) = 9? 
. Therefore, A = z is the MVUB estimator. 


Note: 


First consider the reguarity condition: 


g| ben eld) 7 


Note: 


log p (x|6) log p (x|8) d p (z|6) 
ot = | eas | ae 


Now assuming that we can interchange order of differentiation and integration 


p{ 2s (el0) _ af pede ar | 


00 00 00 : 


So the regularity condition is satisfied whenever this interchange is possible[ footnote]; i.e., when 
derivative is well-defined, fails for uniform density. 

This is simply the Fundamental Theorem of Calculus applied to p (x|0). So long as p (2|6) is 
absolutely continuous with respect to the Lebesgue measure, this is possible. 


Now lets derive the CRLB for a scalar parameter a = g(0), where the pdf is p (x |@). Consider any 
unbiased estimator of a: 


a € (E[a] = a = g(8)) 


Note that this is equivalent to 
[ar (e\a az =o) 


where @ is unbiased. Now differentiate both side 


fee a 
00 00 


or 


— (x |0) 
a 


70 p (z|0) dz = —— 


Now, exmploiting the regularity condition, 
Equation: 


[ @-4) log Pte?) COC 0g(8) 


since 


forex” p (z|0) dz = aEllog p (x|0)| = 0 


Now apply the Cauchy-Schwarz inequality to the integral above: 


(240)'- (J 0) MRE ocanas) 


(220.)" < [ @-a)* velo ae f (Seer) a 


o(a)*is f (a — a)” p (x|0) dz, so 
Equation: 


ors C8) 


- Alogp(z|8) \? 
o|(ee) | 
Now we note that 


E 


(age) _ p| eer el) 


Why? Regularity condition. 


dlogp (x|9)] _ f Plog p (x|9) _ 
Thus, 


of Cows!) 5 (2) da 
00 


=0 


or 


dlog p (x|9) O p (x|9) 
0g oo 


dzx=0 


O° log p (x|) | 
| aan? 219° 


Therefore, 


= ple |0) )| 


O7log p (x|0) 


0 


Thus, [link] becomes 


Jlog p (x|@) Alog p (x|9) 
= 0)dz=E 
|= [7 PEND 5 (a\0) dz 


Note:If 9(9) = 0, then numerator is 1. 


Example:DC Level in White Guassian Noise 
Yn, €4{1,..2,N\:(@, =A + w,) 


where 
Wn ~ AO, a”) 

1 (ae Daa 4") 
(2707) 


#)_ 1% (2, — A)? 
Bel ad a ta eS 


p («|A) = 


n|2 


log p (x| A) 
OA 
log p (x| A) 

| ee eee — 

ecuneale? 

07 log p (z| A) N 


Omsup =o? 


Therefore, the variance of any unbiased estimator satisfies: 
2 


o(a). z= = 


The sample-mean estimator A = ~~ ae Xp, attains this bound and therefore is MVUB. 


Corollary 


When the CRLB is attained 


(8) = Ta 


where 


16) - a 


Omsup 
The quantity I(@) is called Fisher Information that x contains about 0. 


By CRLB Theorem, 


a\ 2 1 
o(8) ~ —E| “ee | 
and 
Plog P(e) ~ 1(6) (6 7 6) 
This yields 


which in turn yields 


So, 


The CRLB is not always attained. 


Example: 
Phase Estimation 


Yn,n € {1,...,N}: (an = Acos(2rfon + vy) + wn) 
The amplitude and frequency are assumed known 


Wn ~ VAG, a”) 


p (a | ~) — a | aa In—A cos(2fon+y)) 
(2n02)? 
01 nA , 
oe 7 ( = dn sin(2mfon + y) — > sin(4r fon + ¢) 


a7] ~ 
< P (zI¢) ( Z ) S> Ln cos(2a fon + yp) — Acos(2rfon + 2y) 


n=l 
d7log p(z|y)]_ A? 
—E = = — 1/2 + 1/2 cos(4mfon + 2y) — cos(4mfon + 2) 


= 
Since I(y) = —B| Salo) |, 


0 


NA? 
I ~ 
Oa 


because Vfp,0 < fo <k: (+ > cos(4r fon) ~ 0) Therefore, 


20? 
—————— 
—~ NA2 


In this case, it can be shown that there does not exist a g such that 


Jlog p (x|¢) 


Bs # I(y) (g(x) — ¢) 


Therefore, an unbiased phase estimator that attains the CRLB does not exist. 
However, a MVUB estimator may still exist--only its variance will be larger than the CRLB. 


Efficiency 


An estimator which is unbiased and attains the CRLB is said to be efficient. 


Example: 
Sample-mean estimator is efficient. 


Example: 

Supposed three unbiased estimators exist for a param 0. 
[missing_resource: ] 

[missing resource: ] 


Example: 
Sinusoidal Frequency Estimation 


Vf0,0 < fo < 1/2: (5n(fo) = Acos(2rfon + ¢)) 
Vn,n S 4h; oe Tee e (tn = 8n(fo) + Wn) 
A and y are known, while fo is unknown. 
y 


Oo 


eS) 
= 
o(fo) ~ A? Ssuue (27n sin(2a fon + y)) 


Suppose _ = 1 (SNR), where N = 10 and y = 0. 
[missing resource: | 


Some 
frequencies are 
easier to 
estimator 
(lower CRLB, 
but not 
necessarily just 
lower bound) 
than others. 


CRLB for Vector Parameter 


0 is unbiased, i.e., 
Vi,i € {1,...,p}: @a » 6,) 
CRLB 


0 (4,) : (I())* 


IV 


where 


d*log p (#8) 
ViANgj: | 1(0),.= —E 
v7 J ( ( dig 00; 00; 
I(0) is the Fisher Information Matrix. 
Cramer-Rao Lower Bound - Vector Parameter 


Assume the pdf p (a |) satisfies the "regularity" condition 


(oe -9 


Then the convariance matrix of any unbiased estimator @ satisfies 
~ (1(8))* >0 


(meaning C, — (I (@))* is p.s.d.) The Fisher Information matrix is 


d*log p (#|8) 
Omsup 


Furthermore, 6 attains the CRLB ( C,= (1 (0))~*) iff 


1(8); ,= B| 


Olog p (x|@) 
——a9 = (9) (9() — 8) 


and 


Example:DC Level in White Guassian Noise 
Yn, mee {1.1 Nia, = A+ w,) 


A is unknown and w,, ~ VY (0, o”) , where o? is unknown. 


!) 


Olog p (a|6) N it 2 
omy 9 a ae 


(Gee a N 


Omsup o? Go 
07 log p («| 6) le 
( @Ado® —*\ yao a 
07 log p (| 0) N he 2 N 
= ee | eee 
( Omsup 20 Dee le 204 
Which leads to 
ne=(" 3 
0 ag 
BND Ge 
Soyein 
6(4) 25 
—\2 _ 204 
2 Sy aS 
o(o4) 2 


Note that the CRLB for A is the same whether or not o? is known. This happens in this case due to 
the diagonal nature of the Fisher Information Matrix. 


In general the Fisher Information Matrix is not diagonal and consequently the CRLBs will depend 
on other unknown parameters. 


Glossary 


idd 
independent and identically distributed 


Maximum Likelihood Estimation 


The maximum likelihood estimator (MLE) is an alternative to the 
minimum variance unbiased estimator (MVUE). For many estimation 
problems, the MVUE does not exist. Moreover, when it does exist, there is 
no systematic procedure for finding it. In constrast, the MLE does not 
necessarily satisfy any optimality criterion, but it can almost always be 
computed, either through exact formulas or numerical techniques. For this 
reason, the MLE is one of the most common estimation procedures used in 
practice. 


The MLE is an important type of estimator for the following reasons: 


1. The MLE implements the likelihood principle. 

2. MLEs are often simple and easy to compute. 

3. MLEs have asymptotic optimality properties (consistency and 
efficiency). 

4. MLEs are invariant under reparameterization. 

5. If an efficient estimator exists, it is the MLE. 

6. In signal detection with unknown parameters (composite hypothesis 
testing), MLEs are used in implementing the generalized likelihood 
ratio test (GLRT). 


This module will discuss these properties in detail, with examples. 


The Likelihood Principle 


Supposed the data X is distributed according to the density or mass 
function p (a|@). The likelihood function for 0 is defined by 


| (O|a) =p (x8) 


At first glance, the likelihood function is nothing new - it is simply a way of 
rewriting the pdf/pmf of X. The difference between the likelihood and the 
pdf or pmf is what is held fixed and what is allowed to vary. When we talk 
about the likelihood, we view the observation z as being fixed, and the 
parameter @ as freely varying. 


Note: It is tempting to view the likelihood function as a probability density 
for 9, and to think of 1 (@|a) as the conditional density of 8 given x. This 
approach to parameter estimation is called fiducial inference, and is not 
accepted by most statisticians. One potential problem, for example, is that 
in many cases ] (@| a) is not integrable ( [ 1 (@|a) d 6 — oo) and thus 
cannot be normalized. A more fundamental problem is that 8 is viewed as 
a fixed quantity, as opposed to random. Thus, it doesn't make sense to talk 
about its density. For the likelihood to be properly thought of as a density, a 
Bayesian approach is required. 


The likelihood principle effectively states that all information we have 
about the unknown parameter @ is contained in the likelihood function. 
PrincipleLikelihood Principle 


The information brought by an observation x about @ is entirely contained 
in the likelihood function p (a|0). Moreover, if x; and x2 are two 
observations depending on the same parameter 9, such that there exists a 
constant c satisfying p (1, |0) = c p (x_|9) for every 9, then they bring 
the same information about @ and must lead to identical estimators. 

In the statement of the likelihood principle, it is not assumed that the two 
observations x; and x2 are generated according to the same model, as long 
as the model is parameterized by @. 


Example: 

Suppose a public health official conducts a survey to estimate 0 < # < 1, 
the percentage of the population eating pizza at least once per week. As a 
result, the official found nine people who had eaten pizza in the last week, 
and three who had not. If no additional information is available regarding 
how the survey was implemented, then there are at least two probability 
models we can adopt. 


1. The official surveyed 12 people, and 9 of them had eaten pizza in the 
last week. In this case, we observe x; = 9, where 


xz, ~ Binomial (12, 6) 


The density for x1 is 
12 
f (0/8) = (Joc — 8) 
L1 


2. Another reasonable model is to assume that the official surveyed 
people until he found 3 non-pizza eaters. In this case, we observe 
2 = 12, where 


Ly ~ NegativeBinomial (3, 1 — 0) 


The density for x2 is 


g (2216) = (7, Jos 0)" 


The likelihoods for these two models are proportional: 
0(0|x1) x £(0|x2) « 69(1 — 6)° 


Therefore, any estimator that adheres to the likelihood principle will 
produce the same estimate for 0, regardless of which of the two data- 
generation models is assumed. 


The likelihood principle is widely accepted among statisticians. In the 
context of parameter estimation, any reasonable estimator should conform 
to the likelihood principle. As we will see, the maximum likelihood 
estimator does. 


Note: While the likelihood principle itself is a fairly reasonable assumption, 
it can also be derived from two somewhat more intuitive assumptions 


known as the sufficiency principle and the conditionality principle. See 
Casella and Berger, Chapter 6. 


The Maximum Likelihood Estimator 


The maximum likelihood estimator 0(a) is defined by 
0 =argmaxl (O| a) 


Intuitively, we are choosing @ to maximize the probability of occurrence of 
the observation z. 


Note: It is possible that multiple parameter values maximize the likelihood 
for a given z. In that case, any of these maximizers can be selected as the 
MLE. It is also possible that the likelihood may be unbounded, in which 
case the MLE does not exist. 


The MLE rule is an implementation of the likelihood principle. If we have 
two observations whose likelihoods are proportional (they differ by a 
constant that does not depend on @), then the value of @ that maximizes one 
likelihood will also maximize the other. In other words, both likelihood 
functions lead to the same inference about 0, as required by the likelihood 
principle. 


Understand that maximum likelihood is a procedure, not an optimality 
criterion. From the definition of the MLE, we have no idea how close it 
comes to the true parameter value relative to other estimators. In constrast, 
the MVUE is defined as the estimator that satisfies a certain optimality 
criterion. However, unlike the MLE, we have no clear produre to follow to 
compute the MVUE. 


Computing the MLE 


If the likelihood function is differentiable, then @ is found by differentiating 
the likelihood (or log-likelihood), equating with zero, and solving: 
Olog 1 (8| ax) 
00 


If multiple solutions exist, then the MLE is the solution that maximizes 
log 1 (8| az), that is, the global maximizer. 


In certain cases, such as pdfs or pmfs with an esponential form, the MLE 
can be easily solved for. That is, 


Olog 1 (O|a) 
00 


can be solved using calculus and standard linear algebra. 


Example: 

DC level in white Guassian noise 

Suppose we observe an unknown amplitude in white Gaussian noise with 
unknown variance: 


Ln —=~At+wWn 


n€ {0,1,...,N —1}, where w, ~ N (0, o”) are independent and 
identically distributed. We would like to estimate 


(2 


by computing the MLE. Differentiating the log-likelihood gives 


Dlog p (w|6) | als 
OA — g 


Olog p (a|@) N Nee 
Ee? a 


Oo? o 


and 


Note: co? is biased! 


As an exercise, try the following problem: 
Exercise: 


Problem: 


Suppose we observe a random sample # = (21...x nv) of Poisson 
measurements with intensity A: Prix; = n] = gre 


n € {0,1,2,...}. Find the MLE for . 


Unfortunately, this approach is only feasible for the most elementary pdfs 
and pmfs. In general, we may have to resort to more advanced numerical 
maximization techniques: 


1. Newton-Raphson iteration 
2. Iteration by the Scoring Method 
3. Expectation-Maximization Algorithm 


All of these are iterative techniques which posit some initial guess at the 
MLE, and then incrementally update that guess. The iteration procedes until 
a local maximum of the likelihood is attained, although in the case of the 
first two methods, such convergence is not guaranteed. The EM algorithm 
has the advantage that the likelihood is always increased at each iteration, 
and so convergence to at least a local maximum is guaranteed (assuming a 
bounded likelihood). For each algorithm, the final estimate is highly 
dependent on the initial guess, and so it is customary to try several different 
starting values. For details on these algorithms, see Kay, Vol. I. 


Asymptotic Properties of the MLE 


Let @ = (@1...2 Nn)" denote an IID sample of size NV, and each sample is 


distributed according to p (a|@). Let 6y denote the MLE based on a 
sample z. 
Asymptotic Properties of MLE 


If the likelihood £(@|a) =p (a|@) satisfies certain "regularity" 


conditions| footnote], then the MLE @y is consistent, and moreover, 6 


converges in probability to 6, where 
6 ~ (0, (1(8)) *) 


where J(@) is the Fisher Information matrix evaluated at the true value of 
0. 

The regularity conditions are essentially the same as those assumed for the 
Cramer-Rao lower bound: the log-likelihood must be twice differentiable, 


and the expected value of the first derivative of the log-likelihood must be 
zero. 


Since the mean of the MLE tends to the true parameter value, we say the 
MLE is asymptotically unbiased. Since the covariance tends to the inverse 
Fisher information matrix, we say the MLE is asymptotically efficient. 


In general, the rate at which the mean-squared error converges to zero is not 
known. It is possible that for small sample sizes, some other estimator may 
have a smaller MSE.The proof of consistency is an application of the weak 
law of large numbers. Derivation of the asymptotic distribution relies on the 
central limit theorem. The theorem is also true in more general settings 
(e.g., dependent samples). See, Kay, Vol. I, Ch. 7 for further discussion. 


The MLE and Efficiency 


In some cases, the MLE is efficient, not just asymptotically efficient. In 
fact, when an efficient estimator exists, it must be the MLE, as described by 
the following result: 


If 8 is an efficient estimator, and the Fisher information matrix (8) is 


positive definite for all 0, then 6 maximizes the likelihood. 


Recall the 6 is efficient (meaning it is unbiased and achieves the Cramer- 
Rao lower bound) if and only if 


omP (29) ~ 1(6) (6 - 6) 


for all 8 and x. Since 0 is assumed to be efficient, this equation holds, and 
in particular it holds when 89 = 6(a). But then the derivative of the log- 

likelihood is zero at 9 = 6(a). Thus, 9 is a critical point of the likelihood. 
Since the Fisher information matrix, which is the negative of the matrix of 


second order derivatives of the log-likelihood, is positive definite, 6 must 
be a maximum of the likelihood. 


An important case where this happens is described in the following 
subsection. 


Optimality of MLE for Linear Statistical Model 


If the observed data x are described by 
a= HO+w 


where H is N x p with full rank, 0 is p x 1, and w ~ -V(O,C), then the 
MLE of @ is 


6=(H’C 1H) ‘H'C1® 


This can be established in two ways. The first is to compute the CRLB for 0 


. It turns out that the condition for equality in the bound is satisfied, and 0 
can be read off from that condition. 


The second way is to maximize the likelihood directly. Equivalently, we 
must minimize 


(a — H0)'C™! (a — H@) 


with respect to @. Since C~! is positive definite, we can write 
C-!=U?TAU = D'D, where D = AtU, where U is an orthogonal 
matrix whose columns are eigenvectors of C’~!, and A is a diagonal matrix 
with positive diagonal entries. Thus, we must minimize 


(Dx — DH0)' (Dx — DH@) 
But this is a linear least squares problem, so the solution is given by the 


pseudoinverse of DH: 
Equation: 


6 = ((DH)" (DH)) (DH)! (Da) 
= (H'™C"!H) ‘H™C 2 
Exercise: 
Problem: 
Consider X1,..., Xn ~ Ns, a1); where s isa p X 1 unknown 


signal, and a? is known. Express the data in the linear model and find 
the MLE Ss for the signal. 


Invariance of MLE 


Suppose we wish to estimate the function w = W(@) and not @ itself. To 
use the maximum likelihood approach for estimating w, we need an 
expression for the likelihood €(w|a) =p (a|w). In other words, we 
would need to be able to parameterize the distribution of the data by w. If 
W is not a one-to-one function, however, this may not be possible. 
Therefore, we define the induced likelihood 


£(w|a) =maxg {0,W(8@) = w}£(O| a) 


The MLE w is defined to be the value of w that maximizes the induced 
likelihood. With this definition, the following invariance principle is 
immediate. 


Let 6 denote the MLE of @. Then ® = w (8) is the MLE of w = W(6). 


The proof follows directly from the definitions of 6 and @. As an exercise, 
work through the logical steps of the proof on your own. 


Example: 
Let @ = (x1...2y)’ where 


x; ~ Poisson (\) 


Given 2, find the MLE of the probability that 2 ~ Poisson (A) exceeds 
the mean A. 


WO)=Prie>A= >> eAX 
n=|A+1| eu 


where |z| = largest integer < z. The MLE of w is 


where J is the MLE of D: 


Be aware that the MLE of a transformed parameter does not necessarily 
satisfy the asymptotic properties discussed earlier. 
Exercise: 


Problem: 


Consider observations 71,...,2, where x; is a p-dimensional vector 
of the form x; = s + w; where s is an unknown signal and w; are 
independent realizations of white Gaussian noise: 


wi ~ NV (0, 07Ipxp) 


Find the maximum likelihood estimate of the energy E = s* s of the 
unknown signal. 


Summary of MLE 


The likelihood principle states that information brought by an observation x 
about @ is entirely contained in the likelihood function p (a@|0). The 
maximum likelihood estimator is one effective implementation of the 
likelihood principle. In some cases, the MLE can be computed exactly, 
using calculus and linear algebra, but at other times iterative numerical 
algorithms are needed. The MLE has several desireable properties: 


e It is consistent and asymptotically efficient (as NV — oo we are doing 
as well as MVUE). 

e When an efficient estimator exists, it is the MLE. 

e The MLE is invariant to reparameterization. 


Bayesian Estimation 


We are interested in estimating 0 given the observation a. Naturally then, 
any estimation strategy will be based on the posterior distribution p (@| a). 
Furthermore, we need a criterion for assessing the quality of potential 
estimators. 


Loss 


The quality of an estimate 6 is measured by a real-valued loss function: 


L (6, 0). For example, squared error or quadratic loss is simply 
“A Paaee & A 
L(@, 6) 2 (6 = 6) (6 = 6) 


Expected Loss 


Posterior Expected Loss: 


B|L(6,6) |2] = [ (6,8) p (6|x) d@ 


Bayes Risk: 
Equation: 


B|L(@,6)| = sr 


I| 
ad 
a 
tS 
BONS 


| 

a 

aE 
&S 


The "best" or optimal estimator given the data a and under a specified loss is 
given by 


0 =argmin E [z (6, @) || 


Example: 
Bayes MSE 


BMSE ( a= ff ( (0-6) p(O|z) dap (w) da 


Since p (a) > 0 for every aw, minimizing the inner integral 
f (@— £[6])? p (6|z) do = B|L(@, 6) | (where B|L(6,6) | is 


the posterior expected loss) for each z, minimizes the overall BMSE. 
Equation: 


as (0-6) v(ixyaa i a ( (0-6) p(0l2)) a 


00 YA) 


-2f (0-6) p (6|z) a6 


Equating this to zero produces 
6 = [or (0|x) d0= El6|a 


The conditional mean (also called posterior mean) of 6 given a! 


Example: 
Vrioe 4 aN SS (Ace) 
W,~ N (0, a) 
prior for unknown parameter A: 
p (a) = U(—Ao, Ao) 


1 : ‘ 
NC) eee Piece) 


if |A] < Ap 
dA 


Minimum Bayes MSE Estimator: 
Equation: 

A = E[Alz] 
= fap (Alz)dA 


-1 y~N 
Ao 4 1 Per AS (en-A)? 4 


—1 
Si DN, (en-A)? 


er dA 


2Ag (2x0?) 2 
Notes 


1. No closed-form estimator 
UN SSah yy awe oy, Ee Ge aio 


3. For smaller Ag, truncated integral produces an A that is a function of x 
woe edule 
2 
4. As N increases, < decreases and posterior p (A| a) becomes tightly 


clustered about = )> zn. This implies A > +> D>, Zn aS N — 00 
(the data "swamps out" the prior) 


Other Common Loss Functions 


Absolute Error Loss 


(Laplace, 1773) 


Equation: 
E|L(8, 6) || = [%|o-4|p (|x) do 


= f°, (6-6) p (6|x) 46+ f° (6-8) p (6|x) a9 
Using integration-by-parts it can be shown that 


[ (0-9) pomas= PO Lies 


—oCo 


(oe) 


[ (0-4) pele)ae= f PO Salads 


where P (8 < y|a) and P (8 > y|a) are a cumulative distributions. So, 


B|L(0,6) |2| [- Pe<vieyay+ [P@>vlayay 


Take the derivative with respect to 6 implies P (9 < 0| x) =P (9 > 8 |x) 


which implies that the optimal 6 under absolute error loss is posterior 
median. 


‘0-1' Loss 
: 0 if O—8@ 
L(6,6) = : Sie, 
(2,8) ere {979} 


[2 (6.8) 2] = 8 [144,912] = Pr[d #6 | 2] 


which is the probability that 6 ~ @ given x. To minimize '0-1' loss we must 
choose 6 to be the value of 6 with the highest posterior probability, which 
implies 6 ~ 0 with the smallest probability. 


p(8|x) 


The optimal estimator 6 under '0-1' loss is the maximum a posteriori 
(MAP) estimator--the value of 8 where p (6|a) is maximized. 


Wiener Filtering and the DFT 


Connecting the Vector Space and Classical Wiener Filters 
Suppose we observe 
Lc=ytrw 


which are all NV x 1 vectors and where w ~ VW (0, oI if Given z we wish 


to estimate y. Think of y as a signal in additive white noise w. a is a noisy 
observation of the signal. 


Taking a Bayesian approach, put a prior on the signal y: 

y~ V(0, Ry) 
which is independent of noise w. The minimum MSE (MMSE) estimator is 

9 = RyRy *@ 
Under the modeling assumptions above 
Equation: 

Ry = E ly(y 5 w)"| 
= EF lyy" | +E [yw | 


= Elyy’| 
= Ryy 


since & [yw | = 0 and since y and w are zero-mean and independent. 
Equation: 


as. = E|xx"| 
= Bl(y+w)(y+w)"| 
= Elyy’|+ Elyw’] + Elwy"| + Elww"| 
= Ry + Row 


since & [ww | — Ryw. Hence 


om 


y = Ryy(Ryy + Ryw)  @ = Hest 
Where H/,,; is the Wiener filter. Recall the frequency domain case 


Syy(f) 
Dad) 7 Dawe t) 


Now let's look at an actual problem scenario. Suppose that we know a priori 
that the signal y is smooth or lowpass. We can incorporate this prior 
knowledge by carefully choosing the prior covariance Ryy. 


Heat) = 


Recall the DFT 
1 al . kn 
Vk,k =0,...N—-1:|%=——) ‘Ye rw) 


or in vector notation 


Vk,k =0,...,.N—1:(%= ly, uz)) 


(4 ian ian 2h ion NUE ) p 
e e eae EC rece 
where uz = Tif (H dehotes Hermitian 


transpose) 


Note: (uz, uz) = ugtu, =—1, VE ALRAL: ((ux, wr) = uz,u; = 0), 
ie., Vk, k =0,...,N —1: ({uz}) is an orthonormal basis. 


The vector uz spans the subspace corresponding to a frequency band 
centered at frequency f;, = ae ("digital" frequency on [0, 1]). If we know 


that y is lowpass, then 
E|(\\ (sue) (I)?] = E[( % II)? 


should be relatively small (compared to E ( I| (y, Wo) I)” ) for high 


frequencies. 


Let 


ox” = B|((\| (yur) I) 


A lowpass model implies 09? > 012 >... > a”, assuming N even, and 
2 


conjugate symmetry implies V7,7 = 1,..., x ; (on_;” = o;”) 
Furthermore, let's model the DFT coefficients as zero-mean and 
independent 


E[|%] = 0 
2. 
= Ok iPt=k 
B|4YY%| = 
%34] eee 


This completely specifies our prior 
y~ MN (0, Fee) 
Ree=UDU= 


where 


0 a1? 
D =, 
0 O ONn_17 
and 
C=(up wi UNn-1) 
Note: 
YyY—Uly 
is the DFT and 
y= UY 


is the inverse DFT. 


With this prior on y the Wiener filter is 
§ = UDU"(UDU" +071) ‘@ 


Since U is a unitary matrix UU" = I and therefore 
Equation: 


UDU¥(U (D+0I)U#) ‘x 
= UDU®U(D +071) ‘U%« 
= UD(D+o0?I) ‘U4e 


te) 
| 


[footnote] Now take the DFT of both sides 


WY =U"G=D(D+07l) & 


where 2 = U"z and is the DFT of x. Both D and D + o7J are diagonal 
SO 


es d 2 
sr «rec ee 4 
dik + 0? OL? +0? 


Hence the Wiener filter is a frequency (DFT) domain filter 


—~ 


Y, = Ay, 2, 


where 2%, is the k*® DFT coefficient of a and the filter response at digital 
frequency aah is 


Assuming 09? > 017 >... > oN and 

9,9 = 1isew + : (on_;” = aj”). The filter's response is a digital 
lowpass filter! 

A Digital Lowpass Filter! 


2 
Oy +o 


0 1 20 N 


If A, B, C are all invertible, compatible matrices, then 
(ABC) + =C"'B1A-1.U+ =U, (UH) =U. 


Summary of Wiener Filter 
Problem: Observe x = y+ w 
Recover/estimate signal y. 


Classical Wiener Filter (continuous-time): 


Syy(w 
H(w) = yy ( ) 
Syy(w) + Srww(w) 
where y(t) and w(t) are stationary processes. 
Vector Space Wiener Filter: 


i= Ryy(Ryy + co ie 


Wiener Filter and DFT: (Ryw = o7 I iO DE H where U is DFT, 
then 4 is a discrete-time filter whose DFT is given by 
Equation: 


N-1 Qa 
H, = Sar h,er" x 


dk 
di.,.+07 


Here, d;,;, plays the same role as Syy(w). 


Kalman Filters 


The Kalman filter is an important generalization of the Wiener filter. Unlike 
Wiener filters, which are designed under the assumption that the signal and 
noise are stationary, the Kalman filter has the ability to adapt itself to non- 
stationary environments. 


The Kalman filter can be viewed as a sequential minimum MSE estimator 
of a signal in additive noise. If the signal and noise are jointly Gaussian, the 
then Kalman filter is optimal in a minimum MSE sense (minimizes 
expected quadratic loss). 


If the signal and/or noise are non-Gaussian, then the Kalman filter is the 
best linear estimator (linear estimator that minimizes MSE among all 
possible linear estimators). 


Dynamical Signal Models 


Recall the simple DC signal estimation problem. 
Equation: 


Vn, n = {0,...,N—1}: (a, = A+ wz) 


Where A is the unknown DC level and w,, is the white Gaussian noise. A 
could represent the voltage of a DC power supply. We know how to find 
several good estimators of A given the measurements {xo,...,2N_1}. 


In practical situations this model may be too simplistic. the load on the 
power supply may charge over time and there will be other variations due to 
temperature and component aging. 


To account for these variations we can employ a more accurate 
measurement model: 
Equation: 


Vn, n = {0,...,N—1}: (an = An + wn) 


where the voltage A,, is the true voltage at time n. 


Now the estimation problem is significantly more complicated since we 
must estimate {Ao,..., Aw—i}. Suppose that the true voltage A,, does not 
vary too rapidly over time. Then successive samples of A,, will not be too 
different, suggesting that the voltage signal displays a high degree of 
correlation. 


This reasoning suggests that it may be reasonable to regard the sequence 
{Ao,..-, Ani}, as a realization of a correlated (not white) random 
process. Adopting a random process model for A,, allows us to pursue a 
Bayesian approach to the estimation problem ((link]). 

Voltage Varying Over Time 


Using the model in [link], it is easy to verify that the maximum likelihood 
and MVUB esitmators are given by 
Equation: 


Ax = In 
Our estimate is simply the noisy measurements! No averaging takes place, 
so there is no noise reduction. 


Let's look at the example again, [link]. 
True Voltage Varying Over Time 


The voltage A, is varying about an average value of 10V. Assume this 
average value is known and write 
Equation: 


A, =10+ yn 


Where y,, is a zero-mean random process. Now a simple model for y, 
which allows us to specify the correlation between samples is the first- 


order Gauss-Markov prcoess model: 
Equation: 


Vn, 4 1235 (Ga = Oa Ee) 


Where un ~ VY (0, Ou") iid (white Gaussian noise process). To initialize 
the process we take yo to be the realization of a Gaussian random variable: 
yovN (0, a. Un is called the driving or excitation noise. The model 
in [link] is called the dynamical or state model. The current output y,, 
depends only on the state of the system at the previous time, or y,_1, and 
the current input u,, ([link]). 


Y1 = Ayo + Uo 
Equation: 


yY2 = ayit Uy 
= a(ayo+uo)+u1 


= ayo + au, + U2 


n 
Yn = a" *tyo F a unk 
k=1 


Equation: 


Elyn) = a” Elyo] + hey a* Blun—al 


0 
Correlation: 
Equation: 
ElymYn] = El(a™™ yo + Dopey a*tm—z) (a"* yo + jy @'Un—1) | 
= Bla ye”) + By DP atm itn] 
Equation: 


on? if m—k=n-l 
E\Um—-kUn—i| = 
[Um Kun] ‘4 otherwise 


If m > n, then 
Equation: 


n 
ElymYn| = a maial a am i y a 
k=1 


If |a| > 1, then it's obvious that the process diverges (variance — oo). 
This is equivalent to having a pole outside the unit circle shown in [link]. 


1 
u ————— y 


So, let's assume |a| < 1 and hence a stable system. Thus as m and n get 
large 


aes +0 


Now let m — n = T. Then for m and n large we have 
Equation: 


fl ieita| Seog Sa 


2. 2 
ao, 
1—a? 


This shows us how correlated the process is: 
ja| + 1 => heavily correlated (or anticorrelated) 
|a| + 0 => weakly correlated 


How can we use this model to help us in our estimation problem? 


The Kalman Filter 


Let's look at a more general formulation of the problem at hand. Suppose 
that we have a vector-valued dynamical equation 
Equation: 


Yn4+1 = Ay, + bun 


Where y,, is p x 1 dimensional, A is p x p, and bis p x 1. The initial 
state vector is Yo ~ V(0, Ry), where Ro is the covariance matrix and 

Un ~ V (0, ai) iid (white Gaussian excitation noise). This reduces to the 
case we just looked at when p = 1. This model could represent a p™ order 
Gauss-Markov process: 

Equation: 


Yn—-1 = Q1Yn + A2Yn-1 +--+ + ApYn—pti + Un 


Define 
Equation: 


Yn—pt1 


Yn—p+2 
yn — 
Yn-1 
Yn 
Then, 
Equation: 
Yn41 = Ayn + bun, 
0 1 0 ts een 0 
0 0 1 «0 O ilautes 
= + + Un 
moe . . 0 
OF 00" caer ache. 10 1 Yn-1 
Q, QQ «++ «++ GAp-1 Ap Yn 1 


Here A is the state transition matrix. Since y, is a linear combination of 
Gaussian vectors: 
Equation: 


Yn = A’yo +S) A® bunt 
k=1 


We know that y,, is also Gaussian distributed with mean and covariance 
R= ava. | Yn ~ WV(, Ry). The covariance can be recursively 
computed from the basic state equation: 

Equation: 


Rnii = AR, A? + 0,,7bb" 


Assume that measurements of the state are available: 
Equation: 


i = CT yn, + wn 


Where wy, ~ VY (0, is") iid independant of {u,,} (white Gaussian 
observation noise). 


For example, if C = (0...01)", then 
Equation: 


Where z,, is the observation, y,, is the signal, and w,, is the noise. Since our 
model for the signal is Gaussian as well as the observation noise, it follows 
that tn ~ W(0,0n”), where o,? = C?R,C + 0,7 (link). 

Block Diagram 


Kalman first posed the problem of estimating the state of y,, from the 
sequence of measurements 


To derive the Kalman filter we will call upon the Gauss-Markov Theorem. 


First note that the conditional distribution of y, given x, is Gaussian: 
Where Ynin is the conditional mean and Pp|,, is the covariance. 


We know that this is the form of the conditional distribution because y,, 
and x, are jointly Gaussian distributed. 


Note: 
where yy is the signal samples yn, . . -, Yn—p+1, En is the 
observations/measurements %p,..-,%n—p+1, aNd Ynjn is the best (minimum 


MSE) estimator of y, given Zp. 


This is all well and good, but we need to know what the conditional mean 
and covariance are explicitly. So the problem is now to find/compute Yp\n 


and Prin. We can take advantage of the recursive state equation to obtain a 
recursive procedure for this calculation. To begin, consider the "predictor" 
Uninet: 


Un eacd a NANG its Pesca) 


Where yr, is the signal samples, {yn, . . ., Uaepicih Zn—1 is the observations 
{@n—1,-++)Ln—p}, and Ynin—1 is the best min MSE estimator of yp given 
Z,—1. Although we don't know what forms Y,),_1 and P,),_1 have, we do 
know two important facts: 


1. The predictor Y,|,1 acts as a sufficient statistic for y,. That is, we can 
replace £,_1 (the data) with Ynjn—1 (the predictor). In other words, all 


the relevant information in x,_; pertaining to y, is summarized by the 
predictor Valids which is, of course, a function of 7,_1. 

2. The predictor Ynjn—1 and the prediction error €njn—1 = Yn — Ynin—1 
are orthogonal (the orthogonality principle of minimum MSE 
estimators => (E Palw=tenbes | = 0) => error is orthogonal to 
estimator). 


Moreover, 
Equation: 


Yn = Ynin—1 zi En|n—1 


Since all quantities are zero-mean, 
En|n—1 ~ N (0, Pin) 


where P,,;,_1 is the covariance of Yn|Ln—1 and "variability" of y,, about the 
predictor YJnjn—1. Therefore, 


Ynli=l 7 N (0, Lin Pigs) 


Where #,, is the covariance of Y,,. Now suppose that we have the predictor 
Yn|n—1 computed and a new measurement is made: 


Equation: 


th. SC yt, 
= Cr (Gnjn—1 ap €n|n—1) + Wn 


Note: §jn—1, €njn—1, aNd wy are all orthogonal. 


We can express all relevant quantities in the matrix equation 


Equation: 


Un I I 0 En|n—1 
Yn|n—1 = 0 I 0 Yn|n—1 
se Cc Cc A Wn 


Now because of the orthogonality, the covariance 
Equation: 


T 
En|n—-1 En|n—-1 Pile 0 0 
E Yn|n—1 Yn|n—1 — 0 Ry, — Pajn—1 0 


iD g Wn, 0) 0 Ow 


Combining this with the matrix [link] shows that 


Equation: 
Yn Yn T Ry Pain-1 R,C 
E Ynin—1 Ynin—1 — Pant 
Ln Ln 
CTR», 8 


P it ae ' en 
Note: (nga) are jointly Gaussian with the covariance in [link] 
and means zero. 


Where 
Equation: 


Pr\n-1 = Rn - n|n—1 


Equation: 


We now have all the quantities necessary to compute our recursive 
estimator using the Gauss-Markov Theorem. 


We will now derive a recursion for conditional distribution of y,, given 
Yn|n—1 (best estimate based on past observations) and x7, (current 


observation). We know that Yn| ( Upasts Xn) is Gaussian (since all 
quantities are jointly Gaussian). Let's denote this conditional distribution by 


Yn | (Orin i) co AC Gis Pri) 


Applying the Gauss-Markov Theorem we find 
Equation: 


which is the best estimator of y,, given Valet and z,,. The inverse of S,, is 


given by 
Equation: 
a c\ /-cr\" 
art, (Pp) @ tm )(- 
0 0 1 1 
where 


Equation: 


ae = on = ct (Rn _ ini) C 
CT Princ + Ca 


Note: 0,2 = C’?R,C + dy? 


Substituting this inverse formula into [link] yields 
Equation: 


Yn|n = Yn|n—1 = Pin Cy, (ty = CT Gnin-1) 


The Gauss-Markov Theorem also gives us an expression for P,,\,: 
Equation: 


and upon substituting [link] for S,, _’ we get 
Equation: 


= i 
Pain — Pajn-1 — Yn Prin-1 CC Prjn-1 
Note that both expressions contain the quantity 
Equation: 


Ky = n|n—1 Cr¥n 


which is the so-called Kalman gain. 


Using the Kalman gain, the Kalman recursions are given by 
Equation: 


Ynin — Yn|n—1 + Ky (an _ OF tna) 
Equation: 


Prin = £n|n—-1 YnKn Ky” 


The recursions are complete except for definitions of Y,),—1 and Prjn—1- 
Equation: 


Yn|n—1 = Elyn|@n—1| 
= El Ayn—1 + buin—1|2n-1] 
= AGn—1)n-1 


Equation: 


Prjn-1 = El (mn _ Inin-1) (Yn — Injnt) | 
= AP»-1)n-1A* te on bb* 


Now we can summarize the Kalman filter: 


1. [link], where Ynin is the best estimate if y,, given observations up to 


time n. 
2. [link] 
3. [link] 
4. [link] 
5. [link] 
6. [link] 


Measurements/observation model: 


Wn ~ N(0;a%5") 
(C = 1). 


Example: 
First-order Gauss-Markov Process 


Yn4+1 = 4Yn + Un 


Where y,+1 is a time-varying voltage, u, ~ V (0, cre and o,, = 0.1. 
(a = 0.99) => highly correlated process. (A = a, b = 1). 


Signal y 
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Kalman Filtering Equations 


1. Prjn—1 = @? Pr—in—1 + On” (q(n) in MATLAB code) 
2. Yn * = Pain—1 + Ow" (g(r) in MATLAB code) 


3. Pain = Prin—1 — Yn *Pnjn—1” (p(n) in MATLAB code) 

4. Ky = Prin-1Yn” (k(m) in MATLAB code) 

5. Inin—1 = @Yn—1|n—1 (py (nm) in MATLAB code) 

6. Inin = Ynin-1 + Kn (@n — Gnin—1) (ey (nm) in MATLAB code) 


Initialization: ey (1) = 0, q(1) = oy” 


0 100 200 300 400 500 600 700 800 9300 1000 
Observations x 
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Introduction to Adaptive Filtering 


The Kalman filter is just one of many adaptive filtering (or estimation) 
algorithms. Despite its elegant derivation and often excellent performance, 
the Kalman filter has two drawbacks: 


1. The derivation and hence performance of the Kalman filter depends on 
the accuracy of the a priori assumptions. The performance can be less 
than impressive if the assumptions are erroneous. 

2. The Kalman filter is fairly computationally demanding, requiring 
O (2 ) operations per sample. This can limit the utility of Kalman 
filters in high rate real time applications. 


As a popular alternative to the Kalman filter, we will investigate the so- 
called least-mean-square (LMS) adaptive filtering algorithm. 


The principle advantages of LMS are 


1. No prior assumptions are made regarding the signal to be estimated. 
2. Computationally, LMS is very efficient, requiring O(P) per sample. 


The price we pay with LMS instead of a Kalman filter is that the rate of 


convergence and adaptation to sudden changes is slower for LMS than for 
the Kalman filter (with correct prior assumptions). 


Adaptive Filtering Applications 


Channel/System Identification 


Channel/System Identification 


noisy 
output 


Noise Cancellation 


Suppression of maternal ECG component in fetal ECG ([link]). 


[missing resource: .png] 


Cancelling 
maternal heartbeat 
in fetal 
electrocardiograph 
y (ECG): position 
of leads. 


noise 


primary signal 
from abdominal 
sensor 


reference signal 
from chest sensor 
(maternal ECG) 


y is an estimate of the maternal ECG signal present in the abdominal signal 
({link]). 


[missing resource: .png] 


Results of fetal 
ECG experiment 
(bandwidth, 3- 
35Hz; sampling 
rate, 256Hz): 
(a)reference input 
(chest lead); 
(b)primary input 
(abdominal lead); 
(c)noise-canceller 
output. 


Channel Equalization 


Channel Equalization 


Adaptive Controller 


Adaptive Controller 
output input 


adaptive 
ee 


Here, the reference signal is the 
desired output. The adaptive controller 
adjusts the controller gains (filter 
weights) to keep them appropriate to 
the system as it changes over time. 


Iterative Minimization 


Most adaptive filtering alogrithms (LMS included) are modifications of 
standard iterative procedures for solving minimization problems in a real- 
time or on-line fashion. Therefore, before deriving the LMS algorithm we 
will look at iterative methods of minimizing error criteria such as MSE. 


Conider the following set-up: 
x, : observation 


Yk : Signal to be estimated 


Linear estimator 


Equation: 
YR = W1 Ly + WoLg_y +... + WyLh_py 
FIR “a 
* Filter ¥4, 


Impulse response of the filter: 


pena Os O81 AO ice St OO p25. 


p?) 


Vector notation 


Equation: 


Where 


Lk 


Lk-1 
t= 
Lk—pt1 

and 

W1 

Ww2 

w = 

bes 
Error signal 
Equation: 

Ck — Yk— Uk 
= va 
= Yr—-—te Ww 

Assumptions 


(xx, Yx) are jointly stationary with zero-mean. 


MSE 


Equation: 


Ele,”| = Bl (y, — tp w) " 
= Elyx”| — 2wl Elaxpy,] + wt E [open |w 
= R,, —2w'R,, + w? Rw 


Where Ryy is the variance of Yk-, Rexx is the covariance matrix of x;, and 
Ryy = E [z.Yx] is the cross-covariance between x; and y; 


Note:The MSE is quadratic in W which implies the MSE surface is 
"bowl" shaped with a unique minimum point ({link]). 


MSE 


Optimum Filter 


Minimize MSE: 
Equation: 

OE |e; 

(Se = 2Ry + 2Ryxw = 7 > (Woy = Ra Ray) 
Ow 
Notice that we can re-write [link] as 
Equation: 
Elx,x," w| = Elrpyr| 
or 
Equation: 
Elz, (y,— 2,7 w)| = Elazex 
= () 


Which shows that the error signal is orthogonal to the input x; (by the 
orthogonality principle of minimum MSE estimator). 


Steepest Descent 


Although we can easily determine wopt by solving the system of equations 
Equation: 


Ryxw = Ryy 


Let's look at an iterative procedure for solving this problem. This will set 
the stage for our adaptive filtering algorithm. 


We want to minimize the MSE. The idea is simple. Starting at some initial 
weight vector wo, iteratively adjust the values to decrease the MSE ((link]). 
In One-Dimension 


MSE 


We want to move wo towards the optimal vector wWop¢. In order to move in 
the correct direction, we must move downhill or in the direction opposite to 
the gradient of the MSE surface at the point wa. Thus, a natural and simple 
adjustment takes the form 

Equation: 


i dE |e,”| 
wi = wo Seay 


Where p is the step size and tells us how far to move in the negative 
gradient direction ([link]). 


MSE 


/ 
7 
Que /, 
cele, vow jf 
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/, 
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Generalizing this idea to an iterative strategy, we get 


Equation: 
1 OE lex? | 
Wk = Wk-1 — o- aw W=Wh_1 
and we can repeatedly update w: wo, w1,..-., we. Hopefully each 


subsequent wy, is closer to Wop. Does the procedure converge? Can we 
adapt it to an on-line, real-time, dynamic situation in which the signals may 
not be stationary? 


LMS Algorithm Analysis 


Objective 
Minimize instantaneous squared error 
Equation: 
e,"(w) = Cr a x,t w)” 
LMS Algorithm 
Equation: 


Wr = We_-1 + MERCK 


Where wz, is the new weight vector, wz_1 is the old weight vector, and 
[4x iS a Small step in the instantaneous error gradient direction. 


Interpretation in Terms of Weight Error Vector 


Define 
Equation: 


Uk = Wk — Wopt 


Where Wopt is the optimal weight vector and 
Equation: 
fs 
Ek = Yk — Lk Wopt 


where €, is the minimum error. The stochastic difference equation is: 
Equation: 


Up = Lvug_1 + UrRER 


Convergence/Stability Analysis 


Show that (tightness) 
Equation: 


limit max {Pr[|] v, ||> Bl} =0 
Bow 


With probability 1, the weight error vector is bounded for all k. 


Chebyshev's inequality is 


Equation: 
B|(\l v II) 
Py[l| ve = B] < ——— 
and 
Equation: 


Pr ve = B] = Sy ((\l Blea] [I)? + o(0)*) 


where (|| E[v,] ||)° is the squared bias. If (|| E[v;] ||)? + o(v,)? is finite 
for all k, then limit Pr[|| vz || > B] = 0 for all k. 
—0o 


Also, 
Equation: 


o(vp)? =r (Elvxvx" |) 


Therefore o(vp)? is finite if the diagonal elements of [y, = E lunve | are 
bounded. 


Convergence in Mean 


|| E[vz] || 0 as & > oo. Take expectation of [link] using smoothing 
property to simplify the calculation. We have convergence in mean if 


1. R,x is positive definite (invertible). 


2 
2.1 = tut hy 


Bounded Variance 


Show that Jy, = EB lv kU | , the weight vector error covariance is bounded 
for all k. 


Note:We could have E[v,] — 0, but o(vz)” — 00; in which case the 
algorithm would not be stable. 


Recall that it is fairly straightforward to show that the diagonal elements of 


the transformed covariance C, = UI,U* tend to zero if pw < a ( Re) ( 


U is the eigenvector matrix of Rx; Ry, = UDU"). The diagonal 
elements of Ci, were denoted by ,;Vi,¢ = {1,...,p}: (¢ = {1,..., p}). 


Note: o(v,)” = tr (I,) = tr (UTC,U) = tr (Ch) = 7? YR 


Thus, to guarantee boundedness of o(v k)? we heed to show that the 
"steady-state" values yz; — (Yi < 00). 


We showed that 
Equation: 


ps (ato,7) 


YE DC gdh) 


where 0,” = E lex? | , A; is the i*® eigenvalue of Ryx ( 


Nix ses. 220 


2 
Re=O FAS U*), anda = “=. 


Equation: 


We found a sufficient condition for yz that guaranteed that the steady-state 
+; 's (and hence o(v,)”) are bounded: 


ewo|bo 


ee er ee 
a ri 


Where )>7_, Ai = tr (Rxx) is the input vector energy. 
With this choice of ~ we have: 


1. convergence in mean 
2. bounded steady-state variance 


This implies 
Equation: 


limit max {Pr[|| v, ||> Bl} =0 
Bow 


In other words, the LMS algorithm is stable about the optimum weight 
vector Wopt- 


Learning Curve 


Recall that 
Equation: 


7 
k = Uk — Lk Wk-1 


and [link]. These imply 
Equation: 

Ck = Ek — Lk VE-1 
where vz = Wk — Wopt. SO the MSE 
Equation: 


Ele,”| 


E|vp-17 open" vp-1| 

= 0,2 4+ E[E|ug_a  epent vp_1|CnénVn,n < kh: (n < k)]] 
= oe" El vp 1 Ryde 1| 

= 0,2 4+ Eltr (Rexve_1vp-1') | 

= o,2+ tr(Roly-1) 


Where (tr (RyxI 4-1) = 1) > (a =— a). So the limiting MSE is 


Equation: 


——— limit Ele,” | 


Since 0 < c < 1 was required for convergence, €4, > 0-7 so that we see 
noisy adaptation leads to an MSE larger than the optimal 
Equation: 


Ele,”| = E| (yx — 24? Wop)’ | 


= O-g 


To quantify the increase in the MSE, define the so-called misadjustment: 
Equation: 


Mo = fexge 


We would of course like to keep M as small as possible. 


Learning Speed and Misadjustment Trade-off 
Fast adaptation and quick convergence require that we take steps as large as 
possible. In other words, learning speed is proportional to yu; larger 4 means 


faster convergence. How does p affect the misadjustment? 


To guarantee convergence/stability we require 


> |colto 


Ear, 


i=l i( Rx) 


Let's assume that in fact w< Foy so that there is no problem with 
i=1 “2 

convergence. This condition implies uw < oT or 

prA; <1Vi,i = {1,...,p}: (@ = {1,...,p}). From here we see that 


Equation: 


1 eee XG iy eee 
c= — ~—p A\i<1 
2 » 1-pri 2 S : 
This misadjustment 
Equation: 
c Tt ee 
M= eS Xi 
l—e pe 


This shows that larger step size pz leads to larger misadjustment. 
Since we still have convergence in mean, this essentially means that with a 


larger step size we "converge" faster but have a larger variance (rattling) 
about Wopt- 


Summary 
small yz implies 


e small misadjustment in steady-state 
e slow adaptation/tracking 


large ys implies 


e large misadjustment in steady-state 
e fast adaptation/tracking 


Example: 


1 
Wopt = 1 


m~(0(6 4) 


if 
Yk = Lk Wopt + Ek 


Ex, ~ V(0,0.01) 
LMS Algorithm 


0 
WE = We_1 + pxperVk, k > 1: (k > 1), where 
Chi hee Ly Wht 
Learning Curve 
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0 
initialization wo = ( ) and 


us = 0.05 


LMS Learning Curve 
[missing_resource: .png] 


{c= O8 


Comparison of Learning Curves 
[missing resource: .png] 


